!! !! This module computes energy bilan on continental surface !! !! @author Marie-Alice Foujols and Jan Polcher !! @Version : $Revision: 1.24 $, $Date: 2009/01/07 13:39:45 $ !! !! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/enerbil.f90,v 1.24 2009/01/07 13:39:45 ssipsl Exp $ !! IPSL (2006) !! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC !! MODULE enerbil ! routines called : restput, restget ! USE ioipsl ! ! modules used : USE constantes USE constantes_veg USE sechiba_io USE parallel ! USE write_field_p, only : WriteFieldI_p IMPLICIT NONE ! public routines : ! enerbil_main ! enerbil_fusion PRIVATE PUBLIC :: enerbil_main, enerbil_fusion,enerbil_clear ! ! variables used inside enerbil module : declaration and initialisation ! LOGICAL, SAVE :: l_first_enerbil=.TRUE. !! Initialisation has to be done one time CHARACTER(LEN=80), SAVE :: var_name !! To store variables names for I/O ! one dimension array allocated, computed and used in enerbil module exclusively REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: psold !! Old surface dry static energy !! saturated specific humudity for old temperature REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsol_sat !! derivative of satured specific humidity at the old temperature REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: pdqsold REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: psnew !! New surface static energy REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsol_sat_new !! New saturated surface air moisture REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: netrad !! Net radiation REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwabs !! LW radiation absorbed by the surface REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwup !! Long-wave up radiation REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwnet !! Net Long-wave radiation REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: fluxsubli !! Energy of sublimation REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsat_air !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tair !! CONTAINS !! !! Main routine for *enerbil* module !! - called only one time for initialisation !! - called every time step !! - called one more time at last time step for writing _restart_ file !! !! Algorithm: !! - call enerbil_begin for initialisation !! - call enerbil_surftemp for psnew and qsol_sat_new !! - call enerbil_flux for tsol_new, netrad, vevapp, fluxlat and fluxsens !! - call enerbil_evapveg for evaporation and transpiration !! !! @call enerbil_begin !! @call enerbil_surftemp !! @call enerbil_flux !! @call enerbil_evapveg !! SUBROUTINE enerbil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, & & index, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, & & qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, & & cimean, ccanopy, emis, soilflx, soilcap, q_cdrag, humrel, fluxsens, fluxlat, & & vevapp, transpir, gpp, vevapnu, vevapwet, vevapsno, t2mdiag, temp_sol, tsol_rad, & & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id ) ! interface description ! input scalar 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 REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write ! input fields INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer 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 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature in Kelvin REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! zonal wind (m/s) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! north-south wind (m/s) 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) :: pb !! Lowest level pressure REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: valpha !! 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) :: emis !! Emissivity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Soil flux REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! This is the cdrag without the wind multiplied REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration in the canopy REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel !! Relative humidity 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) :: vbetaco2 !! Vegetation resistance to CO2 REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: cimean !! mean Ci ! modified fields REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible chaleur flux REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent chaleur flux REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Tsol_rad REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol !! Soil temperature REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! Surface specific humidity REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp !! Assimilation, gC/m**2 total area. REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: t2mdiag !! 2-meter temperature ! ! LOCAL ! REAL(r_std),DIMENSION (kjpindex) :: epot_air_new, qair_new ! ! do initialisation ! IF (l_first_enerbil) THEN IF (long_print) WRITE (numout,*) ' l_first_enerbil : call enerbil_init ' CALL enerbil_init (kjit, ldrestart_read, kjpindex, index, rest_id, qair, temp_sol, temp_sol_new, & & qsurf, tsol_rad, vevapp, fluxsens, fluxlat, gpp, evapot, evapot_corr) CALL enerbil_var_init (kjpindex, temp_air, t2mdiag) RETURN ENDIF ! ! prepares restart file for the next simulation ! IF (ldrestart_write) THEN IF (long_print) WRITE (numout,*) ' we have to complete restart file with ENERBIL variables ' var_name= 'temp_sol' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol, 'scatter', nbp_glo, index_g) var_name= 'qsurf' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, qsurf, 'scatter', nbp_glo, index_g) var_name= 'evapot' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, evapot, 'scatter', nbp_glo, index_g) var_name= 'evapot_corr' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, evapot_corr, 'scatter', nbp_glo, index_g) var_name= 'tsolrad' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, tsol_rad, 'scatter', nbp_glo, index_g) var_name= 'evapora' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, vevapp, 'scatter', nbp_glo, index_g) var_name= 'fluxlat' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, fluxlat, 'scatter', nbp_glo, index_g) var_name= 'fluxsens' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, fluxsens, 'scatter', nbp_glo, index_g) IF ( control%stomate_watchout .OR. control%ok_co2 ) THEN ! The gpp could in principle be recalculated at the beginning of the run. ! However, we would need several variables that are not stored in the restart files. ! var_name= 'gpp' CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, gpp, 'scatter', nbp_glo, index_g) ENDIF RETURN END IF ! ! ! 1. computes some initialisation: psold, qsol_sat and pdqsold ! CALL enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis) ! ! 2. computes psnew, qsol_sat_new, temp_sol_new, qair_new and epot_air_new ! CALL enerbil_surftemp (kjpindex, dtradia, zlev, emis, & & epot_air, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,& & valpha, vbeta1, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, & & qair_new, epot_air_new) ! ! 3. computes lwup, lwnet, tsol_rad, netrad, qsurf, vevapp, evapot, evapot_corr, fluxlat, fluxsubli and fluxsens ! CALL enerbil_flux (kjpindex, dtradia, emis, temp_sol, rau, u, v, q_cdrag, vbeta, valpha, vbeta1, & & 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) ! ! 4. computes in details evaporation and transpiration ! CALL enerbil_evapveg (kjpindex, dtradia, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, cimean, & & ccanopy, rau, u, v, q_cdrag, qair_new, humrel, vevapsno, vevapnu , vevapwet, transpir, gpp, evapot) ! ! 5. diagnose 2-meter temperatures ! CALL enerbil_t2mdiag (kjpindex, temp_air, t2mdiag) IF ( .NOT. almaoutput ) THEN CALL histwrite(hist_id, 'netrad', kjit, netrad, kjpindex, index) CALL histwrite(hist_id, 'evapot', kjit, evapot, kjpindex, index) CALL histwrite(hist_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index) CALL histwrite(hist_id, 'lwdown', kjit, lwabs, kjpindex, index) CALL histwrite(hist_id, 'lwnet', kjit, lwnet, kjpindex, index) IF ( hist2_id > 0 ) THEN CALL histwrite(hist2_id, 'netrad', kjit, netrad, kjpindex, index) CALL histwrite(hist2_id, 'evapot', kjit, evapot, kjpindex, index) CALL histwrite(hist2_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index) CALL histwrite(hist2_id, 'lwdown', kjit, lwabs, kjpindex, index) CALL histwrite(hist2_id, 'lwnet', kjit, lwnet, kjpindex, index) ENDIF ELSE CALL histwrite(hist_id, 'LWnet', kjit, lwnet, kjpindex, index) CALL histwrite(hist_id, 'Qv', kjit, fluxsubli, kjpindex, index) CALL histwrite(hist_id, 'PotEvap', kjit, evapot_corr, kjpindex, index) IF ( hist2_id > 0 ) THEN CALL histwrite(hist2_id, 'LWnet', kjit, lwnet, kjpindex, index) CALL histwrite(hist2_id, 'Qv', kjit, fluxsubli, kjpindex, index) CALL histwrite(hist2_id, 'PotEvap', kjit, evapot_corr, kjpindex, index) ENDIF ENDIF IF (long_print) WRITE (numout,*) ' enerbil_main Done ' END SUBROUTINE enerbil_main !! Algorithm: !! - dynamic allocation for local array !! SUBROUTINE enerbil_init (kjit, ldrestart_read, kjpindex, index, rest_id, qair, temp_sol, temp_sol_new, & & qsurf, tsol_rad, vevapp, fluxsens, fluxlat, gpp, evapot, evapot_corr) ! interface description ! input scalar INTEGER(i_std), INTENT (in) :: kjit !! Time step number LOGICAL ,INTENT (in) :: ldrestart_read !! Logical for _restart_ file to read 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 ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity ! output scalar ! output fields, they need to initialized somehow for the model forcing ORCHIDEE. ! REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol !! Soil temperature REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! near surface specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Tsol_rad REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible chaleur flux REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent chaleur flux REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp !! Assimilation, gC/m**2 total area. REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: evapot !! Soil Potential Evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: evapot_corr !! Soil Potential Evaporation Correction ! local declaration INTEGER(i_std) :: ier ! initialisation IF (l_first_enerbil) THEN l_first_enerbil=.FALSE. ELSE WRITE (numout,*) ' l_first_enerbil false . we stop ' STOP 'enerbil_init' ENDIF ALLOCATE (psold(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in psold allocation. We stop.We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (qsol_sat(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in qsol_sat allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (pdqsold(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in pdqsold allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (psnew(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in psnew allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (qsol_sat_new(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in qsol_sat_new allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (netrad(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in netrad allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (lwabs(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in lwabs allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (lwup(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in lwup allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (lwnet(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in lwnet allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (fluxsubli(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in fluxsubli allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (qsat_air(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in qsat_air allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ALLOCATE (tair(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tair allocation. We stop. We need kjpindex words = ',kjpindex STOP 'enerbil_init' END IF ! open restart input file done by enerbil_init ! and read data from restart input file for ENERBIL process IF (ldrestart_read) THEN IF (long_print) WRITE (numout,*) ' we have to read a restart file for HYDROLOGIC variables' var_name='temp_sol' CALL ioconf_setatt('UNITS', 'K') CALL ioconf_setatt('LONG_NAME','Surface temperature') CALL restget_p (rest_id, var_name, 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 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. ! CALL setvar_p (temp_sol, val_exp,'ENERBIL_TSURF', 280._r_std) ! var_name= 'qsurf' CALL ioconf_setatt('UNITS', 'g/g') CALL ioconf_setatt('LONG_NAME','near surface specific humidity') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., qsurf, "gather", nbp_glo, index_g) IF ( ALL( qsurf(:) .EQ. val_exp ) ) THEN qsurf(:) = qair(:) ENDIF ! var_name= 'evapot' CALL ioconf_setatt('UNITS', 'mm/d') CALL ioconf_setatt('LONG_NAME','Soil Potential Evaporation') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., evapot, "gather", nbp_glo, index_g) ! var_name= 'evapot_corr' CALL ioconf_setatt('UNITS', 'mm/d') CALL ioconf_setatt('LONG_NAME','Corrected Soil Potential Evaporation') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., evapot_corr, "gather", nbp_glo, index_g) ENDIF ! !Config Key = ENERBIL_EVAPOT !Config Desc = Initial Soil Potential Evaporation !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. ! CALL setvar_p (evapot, val_exp, 'ENERBIL_EVAPOT', 0.0_r_std) IF ( ok_var("evapot_corr") ) THEN CALL setvar_p (evapot_corr, val_exp, 'ENERBIL_EVAPOT', 0.0_r_std) ENDIF ! var_name= 'tsolrad' CALL ioconf_setatt('UNITS', 'K') CALL ioconf_setatt('LONG_NAME','Radiative surface temperature') CALL restget_p (rest_id, var_name, 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 ! ! Set the fluxes so that we have something reasonable and not NaN on some machines ! var_name= 'evapora' CALL ioconf_setatt('UNITS', 'Kg/m^2/dt') CALL ioconf_setatt('LONG_NAME','Evaporation') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vevapp, "gather", nbp_glo, index_g) IF ( ALL( vevapp(:) .EQ. val_exp ) ) THEN vevapp(:) = 0.0_r_std ENDIF ! var_name= 'fluxlat' CALL ioconf_setatt('UNITS', 'W/m^2') CALL ioconf_setatt('LONG_NAME','Latent heat flux') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., fluxlat, "gather", nbp_glo, index_g) IF ( ALL( fluxlat(:) .EQ. val_exp ) ) THEN fluxlat(:) = 0.0_r_std ENDIF ! var_name= 'fluxsens' CALL ioconf_setatt('UNITS', 'W/m^2') CALL ioconf_setatt('LONG_NAME','Sensible heat flux') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., fluxsens, "gather", nbp_glo, index_g) IF ( ALL( fluxsens(:) .EQ. val_exp ) ) THEN fluxsens(:) = 0.0_r_std ENDIF ! ! If we are with STOMATE ! IF ( control%stomate_watchout .OR. control%ok_co2 ) THEN ! The gpp could in principle be recalculated at the beginning of the run. ! However, we would need several variables that are not stored in the restart files. var_name= 'gpp' CALL ioconf_setatt('UNITS', 'gC/m**2/time step') CALL ioconf_setatt('LONG_NAME','Gross primary productivity') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., gpp, "gather", nbp_glo, index_g) IF ( ALL( gpp(:,:) .EQ. val_exp ) ) THEN gpp(:,:) = zero ENDIF ENDIF ! initialises temp_sol_new ! saved in thermosoil var_name='temp_sol_new' CALL ioconf_setatt('UNITS', 'K') CALL ioconf_setatt('LONG_NAME','New Surface temperature') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp_sol_new, "gather", nbp_glo, index_g) IF ( ALL( temp_sol_new(:) .EQ. val_exp ) ) THEN temp_sol_new(:) = temp_sol(:) ENDIF ELSE ! initialises temp_sol_new temp_sol_new(:) = temp_sol(:) ENDIF ELSE ! initialises temp_sol_new temp_sol_new(:) = temp_sol(:) ENDIF IF (long_print) WRITE (numout,*) ' enerbil_init done ' END SUBROUTINE enerbil_init SUBROUTINE enerbil_clear () l_first_enerbil=.TRUE. 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) ! ! open restart input file done by enerbil_init ! and read data from restart input file for ENERBIL process END SUBROUTINE enerbil_clear SUBROUTINE enerbil_var_init (kjpindex, temp_air, t2mdiag) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature in Kelvin ! modified fields ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: t2mdiag !! 2-meter temperature CALL enerbil_t2mdiag (kjpindex, temp_air, t2mdiag) IF (long_print) WRITE (numout,*) ' enerbil_var_init done ' END SUBROUTINE enerbil_var_init !! This routines computes psold, qsol_sat, pdqsold and netrad !! SUBROUTINE enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis) ! interface description INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature in Kelvin 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 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: psold !! Old surface dry static energy !! Saturated specific humudity for old temperature REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat !! Derivative of satured specific humidity at the old temperature REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: pdqsold REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: netrad !! Net radiation ! local declaration INTEGER(i_std) :: ji REAL(r_std), DIMENSION(kjpindex) :: dev_qsol REAL(r_std), PARAMETER :: missing = 999998. ! initialisation ! ! 1. computes psold ! psold(:) = temp_sol(:)*cp_air ! ! 2. computes qsol_sat ! 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 (3,'enerbil_begin', & & 'qsol too high ','','') ENDIF ENDDO ENDIF ENDIF ! ! 3. computes pdqsold ! CALL dev_qsatcalc (kjpindex, temp_sol, pb, dev_qsol) DO ji = 1, kjpindex pdqsold(ji) = dev_qsol(ji) * ( pb(ji)**kappa ) / cp_air ENDDO 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 (3,'enerbil_begin', & & 'pdqsold too high ','','') ENDIF ENDDO ENDIF ENDIF ! ! 4. computes netrad and absorbed LW radiation absorbed at the surface ! lwabs(:) = emis(:) * lwdown(:) netrad(:) = lwdown(:) + swnet (:) - (emis(:) * c_stefan * temp_sol(:)**4 + (un - emis(:)) * lwdown(:)) ! IF (long_print) WRITE (numout,*) ' enerbil_begin done ' END SUBROUTINE enerbil_begin !! This routine computes psnew and qsol_sat_new !! !! Computes the energy balance at the surface with an implicit scheme !! that is connected to the Richtmyer and Morton algorithm of the PBL. !! SUBROUTINE enerbil_surftemp (kjpindex, dtradia, zlev, emis, epot_air, & & petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,& & valpha, vbeta1, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, & & qair_new, epot_air_new) ! interface ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT(in) :: dtradia !! Time step in seconds ! input fields 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 !! REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: valpha !! Resistance coefficient REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow 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) :: psnew !! New surface static energy REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat_new !! New saturated surface air moisture REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qair_new !! New air moisture REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: epot_air_new !! New air 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 ! ! ! Help variables ! ! speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) ! zikt = 1/(rau(ji) * speed * q_cdrag(ji)) zikq = 1/(rau(ji) * speed * q_cdrag(ji)) ! ! ! The first step is to compute the fluxes for the old surface conditions ! ! sensfl_old = (petBcoef(ji) - psold(ji)) / (zikt - petAcoef(ji)) larsub_old = chalsu0 * vbeta1(ji) * (peqBcoef(ji) - qsol_sat(ji)) / (zikq - peqAcoef(ji)) lareva_old = chalev0 * (un - vbeta1(ji)) * vbeta(ji) * & & (peqBcoef(ji) - valpha(ji) * qsol_sat(ji)) / (zikq - peqAcoef(ji)) ! ! ! Next the sensitivity terms are computed ! ! netrad_sns = zicp(ji) * quatre * emis(ji) * c_stefan * ((zicp(ji) * psold(ji))**3) sensfl_sns = un / (zikt - petAcoef(ji)) larsub_sns = chalsu0 * vbeta1(ji) * zicp(ji) * pdqsold(ji) / (zikq - peqAcoef(ji)) lareva_sns = chalev0 * (un - vbeta1(ji)) * vbeta(ji) * valpha(ji) * & & zicp(ji) * pdqsold(ji) / (zikq - peqAcoef(ji)) ! ! ! Now we are solving the energy balance ! ! sum_old = netrad(ji) + sensfl_old + larsub_old + lareva_old + soilflx(ji) sum_sns = netrad_sns + sensfl_sns + larsub_sns + lareva_sns dtheta = dtradia * sum_old / (zicp(ji) * soilcap(ji) + dtradia * sum_sns) ! ! psnew(ji) = psold(ji) + dtheta ! qsol_sat_new(ji) = qsol_sat(ji) + zicp(ji) * pdqsold(ji) * dtheta ! temp_sol_new(ji) = psnew(ji) / cp_air ! epot_air_new(ji) = zikt * (sensfl_old - sensfl_sns * dtheta) + psnew(ji) ! fevap = (lareva_old - lareva_sns * dtheta) + (larsub_old - larsub_sns * dtheta) IF ( ABS(fevap) < EPSILON(un) ) THEN qair_new(ji) = qair(ji) ELSE qair_new(ji) = zikq * un / (chalev0 * (un - vbeta1(ji)) * vbeta(ji) * valpha(ji) + & & chalsu0 * vbeta1(ji)) * fevap + qsol_sat_new(ji) ENDIF ! ! ENDDO IF (long_print) WRITE (numout,*) ' enerbil_surftemp done ' END SUBROUTINE enerbil_surftemp !! This routine computes tsol_new, netrad, vevapp, fluxlat, fluxsubli and fluxsens !! SUBROUTINE enerbil_flux (kjpindex, dtradia, emis, temp_sol, rau, u, v, q_cdrag, vbeta, valpha, & & vbeta1, qair, epot_air, psnew, qsurf, fluxsens, fluxlat, fluxsubli, vevapp, temp_sol_new, & & lwdown, swnet, lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT(in) :: dtradia !! Time step in seconds ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature 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 !! REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: valpha !! Resistance coefficient REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: psnew !! New surface static energy REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! Surface specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible chaleur flux REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent chaleur flux REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsubli !! Energy of sublimation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Downward Long-wave radiation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net SW radiation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwup !! Long-wave up radiation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwnet !! Long-wave net radiation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: netrad !! Net radiation REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction ! local declaration INTEGER(i_std) :: ji REAL(r_std),DIMENSION (kjpindex) :: grad_qsat REAL(r_std) :: correction REAL(r_std) :: speed, qc ! initialisation ! ! 1. computes temp_sol_new, netrad, vevapp, fluxlat, fluxsubli, fluxsens ! DO ji=1,kjpindex speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) qc = speed * q_cdrag(ji) 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)) !! !! Add the reflected LW radiation !! lwup(ji) = lwup(ji) + (un - emis(ji)) * lwdown(ji) !! tsol_rad(ji) = (lwup(ji)/ (emis(ji) * c_stefan)) **(1./quatre) !! Need to check the equations !! tsol_rad(ji) = emis(ji) * c_stefan * temp_sol(ji)**4 + lwup(ji) !! !! This is a simple diagnostic which will be used by the GCM to compute the dependence of !! of the surface layer stability on moisture. !! qsurf(ji) = vbeta1(ji) * qsol_sat_new(ji) + & & (un - vbeta1(ji)) * vbeta(ji) * valpha(ji) * qsol_sat_new(ji) qsurf(ji) = MAX(qsurf(ji), qair(ji)) netrad(ji) = lwdown(ji) + swnet(ji) - lwup(ji) vevapp(ji) = dtradia * rau(ji) * qc * vbeta1(ji) * (qsol_sat_new(ji) - qair(ji)) + & & dtradia * rau(ji) * qc * (un - vbeta1(ji)) * vbeta(ji) * & & (valpha(ji) * qsol_sat_new(ji) - qair(ji)) fluxlat(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) *& & (qsol_sat_new(ji) - qair(ji)) + & & chalev0 * rau(ji) * qc * (un - vbeta1(ji)) * vbeta(ji) * & & (valpha(ji) * qsol_sat_new(ji) - qair(ji)) fluxsubli(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) *& & (qsol_sat_new(ji) - qair(ji)) fluxsens(ji) = rau(ji) * qc * (psnew(ji) - epot_air(ji)) lwnet(ji) = lwdown(ji) - lwup(ji) evapot(ji) = MAX(zero, dtradia * rau(ji) * qc * (qsol_sat_new(ji) - qair(ji))) tair(ji) = epot_air(ji) / cp_air ENDDO ! define qsat_air avec subroutine de 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) ! * dtradia !- Penser a sortir evapot en meme temps qu'evapot_corr tdo. DO ji=1,kjpindex speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) qc = speed * q_cdrag(ji) IF ((evapot(ji) .GT. zero) .AND. ((psnew(ji) - epot_air(ji)) .NE. zero )) THEN 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) ) IF (ABS(correction) .GT. min_sechiba) THEN correction = chalev0 * rau(ji) * qc * grad_qsat(ji) * (un - vevapp(ji)/evapot(ji)) / correction ELSE WRITE(numout,*) "Denominateur de la correction de milly nul! Aucune correction appliquee" ENDIF ELSE correction = zero ENDIF correction = MAX (zero, correction) evapot_corr(ji) = evapot(ji) / (un + correction) ENDDO IF (long_print) WRITE (numout,*) ' enerbil_flux done ' END SUBROUTINE enerbil_flux !! This routine computes evaporation and transpiration !! SUBROUTINE enerbil_evapveg (kjpindex, dtradia, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, cimean, & & ccanopy, rau, u, v, q_cdrag, qair, humrel, vevapsno, vevapnu , vevapwet, transpir, gpp, evapot) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT(in) :: dtradia !! Time step in seconds ! input fields 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) :: rau !! Density REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration in the canopy REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel !! Relative humidity 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) :: vbetaco2 !! Vegetation resistance to CO2 REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: cimean !! mean Ci ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp !! Assimilation, gC/m**2 total area. REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception ! local declaration INTEGER(i_std) :: ji, jv REAL(r_std), DIMENSION(kjpindex) :: xx REAL(r_std) :: speed ! initialisation ! ! 1. computes vevapsno, vevapnu ! DO ji=1,kjpindex speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) ! ! 1.1 snow sublimation ! vevapsno(ji) = vbeta1(ji) * dtradia * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji)) ! ! 1.2 bare soil evaporation ! vevapnu(ji) = (un - vbeta1(ji)) * vbeta4(ji) * dtradia * rau(ji) * speed * q_cdrag(ji) & & * (qsol_sat_new(ji) - qair(ji)) END DO ! ! 2. computes transpir and vevapwet ! DO ji = 1, kjpindex ! speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) ! xx(ji) = dtradia * (un-vbeta1(ji)) * (qsol_sat_new(ji)-qair(ji)) * rau(ji) * speed * q_cdrag(ji) ! ENDDO ! DO jv=1,nvm DO ji=1,kjpindex ! ! 2.1 Interception loss ! vevapwet(ji,jv) = xx(ji) * vbeta2(ji,jv) ! ! 2.2 Transpiration ! transpir (ji,jv) = xx(ji) * vbeta3(ji,jv) ! END DO END DO ! ! 3 STOMATE: Assimilation ! IF ( control%ok_co2 ) THEN DO jv = 1, nvm DO ji = 1, kjpindex speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) gpp(ji,jv) = vbetaco2(ji,jv) * dtradia * rau(ji) * speed * q_cdrag(ji) * & (ccanopy(ji) - cimean(ji,jv)) * 12.e-6 ENDDO ENDDO ELSEIF ( control%stomate_watchout ) THEN gpp(:,:) = 0.0 ENDIF IF (long_print) WRITE (numout,*) ' enerbil_evapveg done ' END SUBROUTINE enerbil_evapveg !! Second part of main routine for enerbil module !! - called every time step !! !! Algorithm: !! computes new soil temperature due to ice and snow melt !! SUBROUTINE enerbil_fusion (kjpindex, dtradia, tot_melt, soilcap, temp_sol_new, fusion ) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT(in) :: dtradia !! Time step in seconds ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tot_melt !! Total melt REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity ! modified fields REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! New soil temperature ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fusion !! Fusion ! local declaration INTEGER(i_std) :: ji ! initialisation IF (long_print) WRITE (numout,*) ' enerbil_fusion start ', MINVAL(soilcap), MINLOC(soilcap),& & MAXVAL(soilcap), MAXLOC(soilcap) ! ! 1. computes new soil temperature due to ice and snow melt ! DO ji=1,kjpindex fusion(ji) = tot_melt(ji) * chalfu0 / dtradia temp_sol_new(ji) = temp_sol_new(ji) - ((tot_melt(ji) * chalfu0) / soilcap(ji)) END DO IF (long_print) WRITE (numout,*) ' enerbil_fusion done ' END SUBROUTINE enerbil_fusion !! Diagnose 2 meter air temperature !! SUBROUTINE enerbil_t2mdiag (kjpindex, temp_air, t2mdiag) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature in Kelvin ! modified fields ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: t2mdiag !! 2-meter temperature t2mdiag(:) = temp_air(:) IF (long_print) WRITE (numout,*) ' enerbil_t2mdiag done ' END SUBROUTINE enerbil_t2mdiag END MODULE enerbil