!! !! This module computes hydrologic processes on continental points. !! !! @author Marie-Alice Foujols and Jan Polcher !! @Version : $Revision: 45 $, $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $ !! !< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/hydrolc.f90 $ !< $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $ !< $Author: mmaipsl $ !< $Revision: 45 $ !! IPSL (2006) !! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC !! MODULE hydrolc ! ! ! routines called : restput, restget ! USE ioipsl ! ! modules used : USE constantes USE pft_parameters USE sechiba_io USE grid USE parallel ! USE Write_Field_p IMPLICIT NONE ! public routines : ! hydrol PRIVATE PUBLIC :: hydrolc_main,hydrolc_clear ! ! variables used inside hydrol module : declaration and initialisation ! LOGICAL, SAVE :: l_first_hydrol=.TRUE. !! Initialisation has to be done one time ! LOGICAL, SAVE :: check_waterbal=.FALSE. !! The check the water balance CHARACTER(LEN=80) , SAVE :: var_name !! To store variables names for I/O ! one dimension array allocated, computed, saved and got in hydrol module REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: bqsb !! Hauteur d'eau dans le reservoir profond REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gqsb !! Hauteur d'eau dans le reservoir de surface REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dsg !! Hauteur du reservoir de surface REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dsp !! Hauteur au dessus du reservoir profond REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mean_bqsb !! diagnostique du reservoir profond REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mean_gqsb !! diagnostique du reservoir de surface REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_beg !! Total amount of water at start of time step REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_end !! Total amount of water at end of time step REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_beg !! Total amount of water on vegetation at start of time step REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_end !! Total amount of water on vegetation at end of time step REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_beg !! Total amount of water in the soil at start of time step REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_end !! Total amount of water in the soil at end of time step REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_beg !! Total amount of snow at start of time step REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_end !! Total amount of snow at end of time step REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delsoilmoist !! Change in soil moisture REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delintercept !! Change in interception storage REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delswe !! Change in SWE ! one dimension array allocated, computed and used in hydrol module exclusively REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dss !! Hauteur au dessus du reservoir de surface REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: hdry !! Dry soil heigth in meters REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: precisol !! Eau tombee sur le sol REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: subsnowveg !! Sublimation of snow on vegetation REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snowmelt !! Quantite de neige fondue REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: icemelt !! Quantite de glace fondue REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gdrainage !! Drainage between reservoirs REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vegtot !! Total vegetation ! The last vegetation map which was used to distribute the reservoirs REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: resdist !! Distribution of reservoirs !! profondeur du reservoir contenant le maximum d'eau REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mx_eau_var REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: ruu_ch !! Quantite d'eau maximum REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: runoff !! Ruissellement ! Ajout Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin ! Modifs stabilite REAL(r_std), PARAMETER :: dsg_min = 0.001 CONTAINS !! !! Main routine for *hydrol* 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 hydrolc_snow for snow process (including age of snow) !! - call hydrolc_canop for canopy process !! - call hydrolc_soil for bare soil process !! !! @call hydrolc_snow !! @call hydrolc_canop !! @call hydrolc_soil !! SUBROUTINE hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, & & temp_sol_new, run_off_tot, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,& & qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, tot_melt, transpir, & & precip_rain, precip_snow, returnflow, irrigation, humrel, vegstress, rsol, drysoil_frac, evapot, evapot_corr,& & shumdiag, litterhumdiag, soilcap, 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 INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: returnflow !! Routed water which comes back into the soil REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: irrigation !! Water from irrigation returning to soil moisture REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, ... REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+... REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil capacity REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vevapwet !! Interception loss REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI -> infty) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation Correction ! modified fields REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapnu !! Bare soil evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapsno !! Snow evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow !! Snow mass [Kg/m^2] REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow_age !! Snow age REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio !! Water balance on ice, lakes, .. [Kg/m^2] REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio_age !! Snow age on ice, lakes, ... !! We consider that any water on the ice is snow and we only peforme a water balance to have consistency. !! The water balance is limite to + or - 10^6 so that accumulation is not endless REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Relative humidity REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vegstress !! Veg. moisture stress (only for vegetation growth) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintveg !! Water on vegetation due to interception ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: run_off_tot !! Complete runoff REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: drainage !! Drainage REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (inout):: shumdiag !! relative soil moisture REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rsol !! Resistence to bare soil evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! Fraction of visibly dry soil (between 0 and 1) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_melt !! Total melt ! ! local declaration ! REAL(r_std),DIMENSION (kjpindex) :: soilwet !! A temporary diagnostic of soil wetness REAL(r_std),DIMENSION (kjpindex) :: snowdepth !! Depth of snow layer INTEGER(i_std) :: ji,jv REAL(r_std), DIMENSION(kjpindex) :: histvar !! computations for history files ! ! do initialisation ! IF (l_first_hydrol) THEN IF (long_print) WRITE (numout,*) ' l_first_hydrol : call hydrolc_init ' CALL hydrolc_init (kjit, ldrestart_read, kjpindex, index, rest_id, veget, humrel, vegstress, & & snow, snow_age, snow_nobio, snow_nobio_age, qsintveg) CALL hydrolc_var_init (kjpindex, veget, veget_max, rsol, drysoil_frac, mx_eau_var, ruu_ch, shumdiag, litterhumdiag) ! ! If we check the water balance we first save the total amount of water ! IF (check_waterbal) THEN CALL hydrolc_waterbal(kjpindex, index, .TRUE., dtradia, veget, totfrac_nobio, qsintveg, snow, snow_nobio,& & precip_rain, precip_snow, returnflow, irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno,& & run_off_tot, drainage) ENDIF ! IF (almaoutput) THEN CALL hydrolc_alma(kjpindex, index, .TRUE., qsintveg, snow, snow_nobio, soilwet) ENDIF ! ! shared time step ! IF (long_print) WRITE (numout,*) 'hydrolc pas de temps = ',dtradia 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 HYDROLOGIC variables ' var_name= 'humrel' CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, humrel, 'scatter', nbp_glo, index_g) ! var_name= 'vegstress' CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, vegstress, 'scatter', nbp_glo, index_g) ! var_name= 'snow' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, snow, 'scatter', nbp_glo, index_g) ! var_name= 'snow_age' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, snow_age, 'scatter', nbp_glo, index_g) ! var_name= 'snow_nobio' CALL restput_p(rest_id, var_name, nbp_glo, nnobio, 1, kjit, snow_nobio, 'scatter', nbp_glo, index_g) ! var_name= 'snow_nobio_age' CALL restput_p(rest_id, var_name, nbp_glo, nnobio, 1, kjit, snow_nobio_age, 'scatter', nbp_glo, index_g) ! var_name= 'bqsb' CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, bqsb, 'scatter', nbp_glo, index_g) ! var_name= 'gqsb' CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, gqsb, 'scatter', nbp_glo, index_g) ! var_name= 'dsg' CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, dsg, 'scatter', nbp_glo, index_g) ! var_name= 'dsp' CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, dsp, 'scatter', nbp_glo, index_g) ! var_name= 'dss' CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, dss, 'scatter', nbp_glo, index_g) ! var_name= 'hdry' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, hdry, 'scatter', nbp_glo, index_g) ! var_name= 'qsintveg' CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, qsintveg, 'scatter', nbp_glo, index_g) ! var_name= 'resdist' CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, resdist, 'scatter', nbp_glo, index_g) RETURN ! END IF ! ! computes snow ! CALL hydrolc_snow(kjpindex, dtradia, precip_rain, precip_snow, temp_sol_new, soilcap, & & frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, & & tot_melt, snowdepth) ! ! computes canopy ! ! CALL hydrolc_vegupd(kjpindex, veget, ruu_ch, qsintveg, gqsb, bqsb, dsg, dss,dsp, resdist) ! CALL hydrolc_canop(kjpindex, precip_rain, vevapwet, veget, qsintmax, qsintveg, precisol) ! ! computes hydro_soil ! CALL hydrolc_soil(kjpindex, vevapnu, precisol, returnflow, irrigation, tot_melt, mx_eau_var, veget, ruu_ch, transpir,& & gqsb, bqsb, dsg, dss, rsol, drysoil_frac, hdry, dsp, runoff, run_off_tot, drainage, humrel, vegstress, & & shumdiag, litterhumdiag) ! ! computes horizontal diffusion between the water reservoirs ! IF ( ok_hdiff ) THEN CALL hydrolc_hdiff(kjpindex, dtradia, veget, ruu_ch, gqsb, bqsb, dsg, dss, dsp) ENDIF ! ! If we check the water balance we end with the comparison of total water change and fluxes ! IF (check_waterbal) THEN CALL hydrolc_waterbal(kjpindex, index, .FALSE., dtradia, veget, totfrac_nobio, qsintveg, snow, snow_nobio,& & precip_rain, precip_snow, returnflow, irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno,& & run_off_tot, drainage ) ENDIF ! ! If we use the ALMA standards ! IF (almaoutput) THEN CALL hydrolc_alma(kjpindex, index, .FALSE., qsintveg, snow, snow_nobio, soilwet) ENDIF ! IF ( .NOT. almaoutput ) THEN CALL histwrite(hist_id, 'dss', kjit, dss, kjpindex*nvm, indexveg) CALL histwrite(hist_id, 'bqsb', kjit, mean_bqsb, kjpindex, index) CALL histwrite(hist_id, 'gqsb', kjit, mean_gqsb, kjpindex, index) CALL histwrite(hist_id, 'runoff', kjit, run_off_tot, kjpindex, index) CALL histwrite(hist_id, 'drainage', kjit, drainage, kjpindex, index) CALL histwrite(hist_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg) CALL histwrite(hist_id, 'rain', kjit, precip_rain, kjpindex, index) CALL histwrite(hist_id, 'snowf', kjit, precip_snow, kjpindex, index) CALL histwrite(hist_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg) CALL histwrite(hist_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg) CALL histwrite(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg) histvar(:)=zero DO jv = 1, nvm DO ji = 1, kjpindex IF ( vegtot(ji) .GT. zero ) THEN histvar(ji)=histvar(ji)+veget(ji,jv)/vegtot(ji)*MAX((0.1-dss(ji,jv))*mx_eau_eau, zero) ENDIF ENDDO ENDDO CALL histwrite(hist_id, 'mrsos', kjit, histvar, kjpindex, index) histvar(:)=mean_bqsb(:)+mean_gqsb(:) CALL histwrite(hist_id, 'mrso', kjit, histvar, kjpindex, index) histvar(:)=run_off_tot(:)/one_day CALL histwrite(hist_id, 'mrros', kjit, histvar, kjpindex, index) histvar(:)=(run_off_tot(:)+drainage(:))/one_day CALL histwrite(hist_id, 'mrro', kjit, histvar, kjpindex, index) histvar(:)=(precip_rain(:)-SUM(precisol(:,:),dim=2))/one_day CALL histwrite(hist_id, 'prveg', kjit, histvar, kjpindex, index) ELSE CALL histwrite(hist_id, 'Snowf', kjit, precip_snow, kjpindex, index) CALL histwrite(hist_id, 'Rainf', kjit, precip_rain, kjpindex, index) CALL histwrite(hist_id, 'Qs', kjit, run_off_tot, kjpindex, index) CALL histwrite(hist_id, 'Qsb', kjit, drainage, kjpindex, index) CALL histwrite(hist_id, 'Qsm', kjit, tot_melt, kjpindex, index) CALL histwrite(hist_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index) CALL histwrite(hist_id, 'DelSWE', kjit, delswe, kjpindex, index) CALL histwrite(hist_id, 'DelIntercept', kjit, delintercept, kjpindex, index) ! CALL histwrite(hist_id, 'SoilMoist', kjit, tot_watsoil_end, kjpindex, index) CALL histwrite(hist_id, 'SoilWet', kjit, soilwet, kjpindex, index) CALL histwrite(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg) ! CALL histwrite(hist_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index) CALL histwrite(hist_id, 'SubSnow', kjit, vevapsno, kjpindex, index) ! CALL histwrite(hist_id, 'SnowDepth', kjit, snowdepth, kjpindex, index) ! ENDIF IF ( hist2_id > 0 ) THEN IF ( .NOT. almaoutput ) THEN CALL histwrite(hist2_id, 'dss', kjit, dss, kjpindex*nvm, indexveg) CALL histwrite(hist2_id, 'bqsb', kjit, mean_bqsb, kjpindex, index) CALL histwrite(hist2_id, 'gqsb', kjit, mean_gqsb, kjpindex, index) CALL histwrite(hist2_id, 'runoff', kjit, run_off_tot, kjpindex, index) CALL histwrite(hist2_id, 'drainage', kjit, drainage, kjpindex, index) CALL histwrite(hist2_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg) CALL histwrite(hist2_id, 'rain', kjit, precip_rain, kjpindex, index) CALL histwrite(hist2_id, 'snowf', kjit, precip_snow, kjpindex, index) CALL histwrite(hist2_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg) CALL histwrite(hist2_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg) CALL histwrite(hist2_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg) histvar(:)=zero DO jv = 1, nvm DO ji = 1, kjpindex IF ( vegtot(ji) .GT. zero ) THEN histvar(ji)=histvar(ji)+veget(ji,jv)/vegtot(ji)*MAX((0.1-dss(ji,jv))*mx_eau_eau, zero) ENDIF ENDDO ENDDO CALL histwrite(hist2_id, 'mrsos', kjit, histvar, kjpindex, index) histvar(:)=(run_off_tot(:)+drainage(:))/one_day CALL histwrite(hist2_id, 'mrro', kjit, histvar, kjpindex, index) ELSE CALL histwrite(hist2_id, 'Snowf', kjit, precip_snow, kjpindex, index) CALL histwrite(hist2_id, 'Rainf', kjit, precip_rain, kjpindex, index) CALL histwrite(hist2_id, 'Qs', kjit, run_off_tot, kjpindex, index) CALL histwrite(hist2_id, 'Qsb', kjit, drainage, kjpindex, index) CALL histwrite(hist2_id, 'Qsm', kjit, tot_melt, kjpindex, index) CALL histwrite(hist2_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index) CALL histwrite(hist2_id, 'DelSWE', kjit, delswe, kjpindex, index) CALL histwrite(hist2_id, 'DelIntercept', kjit, delintercept, kjpindex, index) ! CALL histwrite(hist2_id, 'SoilMoist', kjit, tot_watsoil_end, kjpindex, index) CALL histwrite(hist2_id, 'SoilWet', kjit, soilwet, kjpindex, index) ! CALL histwrite(hist2_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index) CALL histwrite(hist2_id, 'SubSnow', kjit, vevapsno, kjpindex, index) ! CALL histwrite(hist2_id, 'SnowDepth', kjit, snowdepth, kjpindex, index) ! ENDIF ENDIF ! IF (long_print) WRITE (numout,*) ' hydrolc_main Done ' END SUBROUTINE hydrolc_main !! Algorithm: !! - dynamic allocation for local array !! - _restart_ file reading for HYDROLOGIC variables !! SUBROUTINE hydrolc_init(kjit, ldrestart_read, kjpindex, index, rest_id, veget, humrel, vegstress, & & snow, snow_age, snow_nobio, snow_nobio_age, qsintveg) ! 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,nvm), INTENT (in) :: veget !! Carte de vegetation ! output fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Stress hydrique, relative humidity REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vegstress !! Veg. moisture stress (only for vegetation growth) REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow !! Snow mass [Kg/m^2] REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow_age !! Snow age REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio !! Snow on ice, lakes, ... REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio_age !! Snow age on ice, lakes, ... REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintveg !! Water on vegetation due to interception ! local declaration INTEGER(i_std) :: ier INTEGER(i_std) :: ji,jv,ik REAL(r_std), DIMENSION (kjpindex,nvm) :: zdsp, tmpdss REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: dsp_g REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: zdsp_g REAL(r_std), DIMENSION(kjpindex) :: a_subgrd ! initialisation IF (l_first_hydrol) THEN l_first_hydrol=.FALSE. ELSE WRITE (numout,*) ' l_first_hydrol false . we stop ' STOP 'hydrolc_init' ENDIF ! >> DS : july 2011 add initialisation of sneige and throughfall_by_pft sneige = snowcri/mille throughfall_by_pft = throughfall_by_pft / 100. ! make dynamic allocation with good dimension ! one dimension array allocation with possible restart value ALLOCATE (bqsb(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in bqsb allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF bqsb(:,:) = zero ALLOCATE (gqsb(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in gqsb allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF gqsb(:,:) = zero ALLOCATE (dsg(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in dsg allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF dsg(:,:) = zero ALLOCATE (dsp(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in dsp allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF dsp(:,:) = zero ! one dimension array allocation ALLOCATE (mean_bqsb(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mean_bqsb allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF mean_bqsb(:) = zero ALLOCATE (mean_gqsb(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mean_gqsb allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF mean_gqsb(:) = zero ALLOCATE (dss(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in dss allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF dss(:,:) = zero ALLOCATE (hdry(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in hdry allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF hdry(:) = zero ALLOCATE (precisol(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF precisol(:,:) = zero ALLOCATE (gdrainage(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF gdrainage(:,:) = zero ALLOCATE (subsnowveg(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in subsnowveg allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF subsnowveg(:) = zero ALLOCATE (subsnownobio(kjpindex,nnobio),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in subsnownobio allocation. We stop. We need kjpindex*nnobio words = ', & kjpindex*nnobio STOP 'hydrolc_init' END IF subsnownobio(:,:) = zero ALLOCATE (snowmelt(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in snowmelt allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF snowmelt(:) = zero ALLOCATE (icemelt(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in icemelt allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF icemelt(:) = zero ALLOCATE (mx_eau_var(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mx_eau_var allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF mx_eau_var(:) = zero ALLOCATE (ruu_ch(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in ruu_ch allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF ruu_ch(:) = zero ALLOCATE (vegtot(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in vegtot allocation. We stop. We need kjpindex words = ',kjpindex*nvm STOP 'hydrolc_init' END IF vegtot(:) = zero ALLOCATE (resdist(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in resdist allocation. We stop. We need kjpindex words = ',kjpindex*nvm STOP 'hydrolc_init' END IF resdist(:,:) = zero ALLOCATE (runoff(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in runoff allocation. We stop. We need kjpindex words = ',kjpindex*nvm STOP 'hydrolc_init' END IF runoff(:,:) = zero ! ! If we check the water balance we need two more variables ! IF ( check_waterbal ) THEN ALLOCATE (tot_water_beg(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tot_water_beg allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF ALLOCATE (tot_water_end(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tot_water_end allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF ENDIF ! ! If we use the almaoutputs we need four more variables ! IF ( almaoutput ) THEN ALLOCATE (tot_watveg_beg(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tot_watveg_beg allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF ALLOCATE (tot_watveg_end(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tot_watveg_end allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF ALLOCATE (tot_watsoil_beg(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tot_watsoil_beg allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF ALLOCATE (tot_watsoil_end(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tot_watsoil_end allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF ALLOCATE (delsoilmoist(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in delsoilmoist allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF ALLOCATE (delintercept(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in delintercept. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' END IF ALLOCATE (delswe(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in delswe. We stop. We need kjpindex words = ',kjpindex STOP 'hydrolc_init' ENDIF ALLOCATE (snow_beg(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in snow_beg allocation. We stop. We need kjpindex words =',kjpindex STOP 'hydrolc_init' END IF ALLOCATE (snow_end(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in snow_end allocation. We stop. We need kjpindex words =',kjpindex STOP 'hydrolc_init' END IF ENDIF ! open restart input file done by sechiba_init ! and read data from restart input file for HYDROLOGIC process IF (ldrestart_read) THEN IF (long_print) WRITE (numout,*) ' we have to read a restart file for HYDROLOGIC variables' var_name= 'snow' CALL ioconf_setatt('UNITS', 'kg/m^2') CALL ioconf_setatt('LONG_NAME','Snow mass') CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow, "gather", nbp_glo, index_g) !!$ ! correction for old restart !!$ DO ik=1, kjpindex !!$ if(snow(ik).gt.maxmass_glacier) snow(ik)=maxmass_glacier !!$ ENDDO ! var_name= 'snow_age' CALL ioconf_setatt('UNITS', 'd') CALL ioconf_setatt('LONG_NAME','Snow age') CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow_age, "gather", nbp_glo, index_g) ! var_name= 'snow_nobio' CALL ioconf_setatt('UNITS', 'kg/m^2') CALL ioconf_setatt('LONG_NAME','Snow on other surface types') CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio, "gather", nbp_glo, index_g) !!$ ! correction for old restart !!$ DO ik=1, kjpindex !!$ if(snow_nobio(ik,iice).gt.maxmass_glacier) snow_nobio(ik,iice)=maxmass_glacier !!$ ENDDO ! var_name= 'snow_nobio_age' CALL ioconf_setatt('UNITS', 'd') CALL ioconf_setatt('LONG_NAME','Snow age on other surface types') CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio_age, "gather", nbp_glo, index_g) ! var_name= 'humrel' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Soil moisture stress') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., humrel, "gather", nbp_glo, index_g) ! var_name= 'vegstress' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Vegetation growth moisture stress') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., vegstress, "gather", nbp_glo, index_g) ! var_name= 'bqsb' CALL ioconf_setatt('UNITS', 'kg/m^2') CALL ioconf_setatt('LONG_NAME','Deep soil moisture') CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., bqsb, "gather", nbp_glo, index_g) ! var_name= 'gqsb' CALL ioconf_setatt('UNITS', 'kg/m^2') CALL ioconf_setatt('LONG_NAME','Surface soil moisture') CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., gqsb, "gather", nbp_glo, index_g) ! var_name= 'dsg' CALL ioconf_setatt('UNITS', 'm') CALL ioconf_setatt('LONG_NAME','Depth of upper reservoir') CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dsg, "gather", nbp_glo, index_g) ! var_name= 'dsp' CALL ioconf_setatt('UNITS', 'm') CALL ioconf_setatt('LONG_NAME','Depth to lower reservoir') CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dsp, "gather", nbp_glo, index_g) ! var_name= 'dss' CALL ioconf_setatt('UNITS', 'm') CALL ioconf_setatt('LONG_NAME','Hauteur au dessus du reservoir de surface') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dss, "gather", nbp_glo, index_g) ENDIF ! var_name= 'hdry' CALL ioconf_setatt('UNITS', 'm') CALL ioconf_setatt('LONG_NAME','Dry soil heigth in meters') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., hdry, "gather", nbp_glo, index_g) ENDIF ! var_name= 'qsintveg' CALL ioconf_setatt('UNITS', 'kg/m^2') CALL ioconf_setatt('LONG_NAME','Intercepted moisture') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., qsintveg, "gather", nbp_glo, index_g) ! var_name= 'resdist' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Distribution of reservoirs') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., resdist, "gather", nbp_glo, index_g) ! ! get restart values if non were found in the restart file ! !Config Key = HYDROL_SNOW !Config Desc = Initial snow mass if not found in restart !Config Def = 0.0 !Config Help = The initial value of snow mass 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 (snow, val_exp, 'HYDROL_SNOW', zero) ! !Config Key = HYDROL_SNOWAGE !Config Desc = Initial snow age if not found in restart !Config Def = 0.0 !Config Help = The initial value of snow age 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 (snow_age, val_exp, 'HYDROL_SNOWAGE', zero) ! !Config Key = HYDROL_SNOW_NOBIO !Config Desc = Initial snow amount on ice, lakes, etc. if not found in restart !Config Def = 0.0 !Config Help = The initial value of snow 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 (snow_nobio, val_exp, 'HYDROL_SNOW_NOBIO', zero) ! !Config Key = HYDROL_SNOW_NOBIO_AGE !Config Desc = Initial snow age on ice, lakes, etc. if not found in restart !Config Def = 0.0 !Config Help = The initial value of snow age 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 (snow_nobio_age, val_exp, 'HYDROL_SNOW_NOBIO_AGE', zero) ! !Config Key = HYDROL_HUMR !Config Desc = Initial soil moisture stress if not found in restart !Config Def = 1.0 !Config Help = The initial value of soil moisture stress 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 (humrel, val_exp,'HYDROL_HUMR', un) CALL setvar_p (vegstress, val_exp,'HYDROL_HUMR', un) ! !Config Key = HYDROL_BQSB !Config Desc = Initial restart deep soil moisture if not found in restart !Config Def = DEF !Config Help = The initial value of deep soil moisture 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. Default behaviour is a saturated soil. ! CALL setvar_p (bqsb, val_exp, 'HYDROL_BQSB', mx_eau_eau*dpu_cste) ! !Config Key = HYDROL_GQSB !Config Desc = Initial upper soil moisture if not found in restart !Config Def = 0.0 !Config Help = The initial value of upper soil moisture 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 (gqsb, val_exp, 'HYDROL_GQSB', zero) ! !Config Key = HYDROL_DSG !Config Desc = Initial upper reservoir depth if not found in restart !Config Def = 0.0 !Config Help = The initial value of upper reservoir depth 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 (dsg, val_exp, 'HYDROL_DSG', zero) ! set inital value for dsp if needed ! !Config Key = HYDROL_DSP !Config Desc = Initial dry soil above upper reservoir if not found in restart !Config Def = DEF !Config Help = The initial value of dry soil above upper reservoir 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. The default behaviour !Config is to compute it from the variables above. Should be OK most of !Config the time. ! zdsp(:,:) = dpu_cste - bqsb(:,:) / mx_eau_eau dsp(1,1) = val_exp CALL getin_p('HYDROL_DSP', dsp(1,1)) IF (dsp(1,1) == val_exp) THEN dsp(:,:) = zdsp(:,:) ELSE IF (is_root_prc) & ALLOCATE(zdsp_g(nbp_glo,nvm),dsp_g(nbp_glo,nvm)) CALL gather(zdsp,zdsp_g) IF (is_root_prc) & CALL setvar (dsp_g, val_exp, 'HYDROL_DSP', zdsp_g) CALL scatter(dsp_g,dsp) IF (is_root_prc) & DEALLOCATE(zdsp_g, dsp_g) ENDIF ! !Config Key = HYDROL_QSV !Config Desc = Initial water on canopy if not found in restart !Config Def = 0.0 !Config Help = The initial value of moisture on canopy 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 (qsintveg, val_exp, 'HYDROL_QSV', zero) ! tmpdss = dsg - gqsb / mx_eau_eau IF ( ok_var("dss") ) THEN CALL setvar_p (dss, val_exp, 'NO_KEYWORD', tmpdss) ELSE dss(:,:)=tmpdss(:,:) ENDIF IF ( ok_var("hdry") ) THEN IF (MINVAL(hdry) .EQ. MAXVAL(hdry) .AND. MAXVAL(hdry) .EQ. val_exp) THEN a_subgrd(:) = zero DO ji=1,kjpindex IF ( gqsb(ji,1)+bqsb(ji,1) .GT. zero ) THEN ! IF (.NOT. (dsg(ji,1).EQ. zero .OR. gqsb(ji,1).EQ.zero)) THEN ! Ajouts Nathalie - Fred - le 28 Mars 2006 a_subgrd(ji)=MIN(MAX(dsg(ji,1)-dss(ji,1),zero)/dsg_min,un) ! ENDIF ENDIF ! ENDDO ! Correction Nathalie - le 28 Mars 2006 - re-ecriture drysoil_frac/hdry - Fred Hourdin ! revu 22 novembre 2007 hdry(:) = a_subgrd(:)*dss(:,1) + (1.-a_subgrd(:))*dsp(:,1) ENDIF ELSE a_subgrd(:) = zero DO ji=1,kjpindex IF ( gqsb(ji,1)+bqsb(ji,1) .GT. zero ) THEN ! IF (.NOT. (dsg(ji,1).EQ. zero .OR. gqsb(ji,1).EQ.zero)) THEN ! Ajouts Nathalie - Fred - le 28 Mars 2006 a_subgrd(ji)=MIN(MAX(dsg(ji,1)-dss(ji,1),zero)/dsg_min,un) ! ENDIF ENDIF ! ENDDO ! Correction Nathalie - le 28 Mars 2006 - re-ecriture drysoil_frac/hdry - Fred Hourdin ! revu 22 novembre 2007 hdry(:) = a_subgrd(:)*dss(:,1) + (un-a_subgrd(:))*dsp(:,1) ENDIF ! ! There is no need to configure the initialisation of resdist. If not available it is the vegetation map ! IF ( MINVAL(resdist) .EQ. MAXVAL(resdist) .AND. MINVAL(resdist) .EQ. val_exp) THEN resdist = veget ENDIF ! ! Remember that it is only frac_nobio + SUM(veget(,:)) that is equal to 1. Thus we need vegtot ! DO ji = 1, kjpindex vegtot(ji) = SUM(veget(ji,:)) ENDDO ! ENDIF ! ! Where vegetation fraction is zero, set water to that of bare soil. ! This does not create any additional water. ! DO jv = 2, nvm DO ji = 1, kjpindex IF ( veget(ji,jv) .LT. EPSILON(un) ) THEN gqsb(ji,jv) = gqsb(ji,1) bqsb(ji,jv) = bqsb(ji,1) dsg(ji,jv) = dsg(ji,1) dss(ji,jv) = dss(ji,1) dsp(ji,jv) = dsp(ji,1) ENDIF ENDDO ENDDO ! DO ik=1, kjpindex if(snow(ik).gt.maxmass_glacier) then WRITE(numout,*)' il faut diminuer le stock de neige car snow > maxmass_glacier dans restart' snow(ik)=maxmass_glacier endif if(snow_nobio(ik,iice).gt.maxmass_glacier) then WRITE(numout,*)' il faut diminuer le stock de neige car snow_nobio > maxmass_glacier dans restart' snow_nobio(ik,iice)=maxmass_glacier endif ENDDO ! IF (long_print) WRITE (numout,*) ' hydrolc_init done ' ! END SUBROUTINE hydrolc_init ! !------------------------------------- ! SUBROUTINE hydrolc_clear() l_first_hydrol=.TRUE. IF (ALLOCATED (bqsb)) DEALLOCATE (bqsb) IF (ALLOCATED (gqsb)) DEALLOCATE (gqsb) IF (ALLOCATED (dsg)) DEALLOCATE (dsg) IF (ALLOCATED (dsp)) DEALLOCATE (dsp) IF (ALLOCATED (mean_bqsb)) DEALLOCATE (mean_bqsb) IF (ALLOCATED (mean_gqsb)) DEALLOCATE (mean_gqsb) IF (ALLOCATED (dss)) DEALLOCATE (dss) IF (ALLOCATED (hdry)) DEALLOCATE (hdry) IF (ALLOCATED (precisol)) DEALLOCATE (precisol) IF (ALLOCATED (gdrainage)) DEALLOCATE (gdrainage) IF (ALLOCATED (subsnowveg)) DEALLOCATE (subsnowveg) IF (ALLOCATED (subsnownobio)) DEALLOCATE (subsnownobio) IF (ALLOCATED (snowmelt)) DEALLOCATE (snowmelt) IF (ALLOCATED (icemelt)) DEALLOCATE (icemelt) IF (ALLOCATED (mx_eau_var)) DEALLOCATE (mx_eau_var) IF (ALLOCATED (ruu_ch)) DEALLOCATE (ruu_ch) IF (ALLOCATED (vegtot)) DEALLOCATE (vegtot) IF (ALLOCATED (resdist)) DEALLOCATE (resdist) IF (ALLOCATED (runoff)) DEALLOCATE (runoff) IF (ALLOCATED (tot_water_beg)) DEALLOCATE (tot_water_beg) IF (ALLOCATED (tot_water_end)) DEALLOCATE (tot_water_end) IF (ALLOCATED (tot_watveg_beg)) DEALLOCATE (tot_watveg_beg) IF (ALLOCATED (tot_watveg_end)) DEALLOCATE (tot_watveg_end) IF (ALLOCATED (tot_watsoil_beg)) DEALLOCATE (tot_watsoil_beg) IF (ALLOCATED (tot_watsoil_end)) DEALLOCATE (tot_watsoil_end) IF (ALLOCATED (delsoilmoist)) DEALLOCATE (delsoilmoist) IF (ALLOCATED (delintercept)) DEALLOCATE (delintercept) IF (ALLOCATED (snow_beg)) DEALLOCATE (snow_beg) IF (ALLOCATED (snow_end)) DEALLOCATE (snow_end) IF (ALLOCATED (delswe)) DEALLOCATE (delswe) ! END SUBROUTINE hydrolc_clear !! This routine initializes HYDROLOGIC variables !! - mx_eau_var !! - ruu_ch !! SUBROUTINE hydrolc_var_init (kjpindex, veget, veget_max, rsol, drysoil_frac, mx_eau_var, ruu_ch, shumdiag, litterhumdiag) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rsol !! Resistance to bare soil evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! Fraction of visible dry soil !! Profondeur du reservoir contenant le maximum d'eau REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: mx_eau_var !! REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ruu_ch !! Quantite d'eau maximum REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out):: shumdiag !! relative soil moisture REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity ! local declaration INTEGER(i_std) :: ji,jv, jd REAL(r_std), DIMENSION(kjpindex) :: mean_dsg REAL(r_std) :: gtr, btr REAL(r_std), DIMENSION(nbdl+1) :: tmp_dl ! ! initialisation tmp_dl(1) = 0 tmp_dl(2:nbdl+1) = diaglev(1:nbdl) ! mx_eau_var(:) = zero ! DO ji = 1,kjpindex DO jv = 1,nvm mx_eau_var(ji) = mx_eau_var(ji) + veget(ji,jv)*wmax_veg(jv)*dpu_cste END DO IF (vegtot(ji) .GT. zero) THEN mx_eau_var(ji) = mx_eau_var(ji)/vegtot(ji) ELSE mx_eau_var(ji) = mx_eau_eau*dpu_cste ENDIF ruu_ch(ji) = mx_eau_var(ji) / dpu_cste END DO ! ! ! could be done with SUM instruction but this kills vectorization mean_bqsb(:) = zero mean_gqsb(:) = zero mean_dsg(:) = zero DO jv = 1, nvm DO ji = 1, kjpindex mean_bqsb(ji) = mean_bqsb(ji) + resdist(ji,jv)*bqsb(ji,jv) mean_gqsb(ji) = mean_gqsb(ji) + resdist(ji,jv)*gqsb(ji,jv) mean_dsg(ji) = mean_dsg(ji) + resdist(ji,jv)*dsg(ji,jv) ENDDO ENDDO mean_dsg(:) = MAX( mean_dsg(:), mean_gqsb(:)/ruu_ch(:) ) DO ji = 1, kjpindex IF (vegtot(ji) .GT. zero) THEN mean_bqsb(ji) = mean_bqsb(ji)/vegtot(ji) mean_gqsb(ji) = mean_gqsb(ji)/vegtot(ji) mean_dsg(ji) = mean_dsg(ji)/vegtot(ji) ENDIF ENDDO DO jd = 1,nbdl ! DO ji = 1,kjpindex !!$ ! !!$ DO jd = 1,nbdl IF ( tmp_dl(jd+1) .LT. mean_dsg(ji)) THEN shumdiag(ji,jd) = mean_gqsb(ji)/mx_eau_var(ji) ELSE IF ( tmp_dl(jd) .LT. mean_dsg(ji)) THEN gtr = (mean_dsg(ji)-tmp_dl(jd))/(tmp_dl(jd+1)-tmp_dl(jd)) btr = 1 - gtr shumdiag(ji,jd) = gtr*mean_gqsb(ji)/mx_eau_var(ji) + & & btr*mean_bqsb(ji)/mx_eau_var(ji) ELSE shumdiag(ji,jd) = mean_bqsb(ji)/mx_eau_var(ji) ENDIF ENDIF shumdiag(ji,jd) = MAX(MIN(shumdiag(ji,jd), un), zero) ENDDO ENDDO ! The fraction of soil which is visibly dry (dry when dss = 0.1 m) drysoil_frac(:) = MIN(MAX(dss(:,1),zero)*10._r_std, un) ! ! Compute the resistance to bare soil evaporation ! rsol(:) = -un DO ji = 1, kjpindex IF (veget(ji,1) .GE. min_sechiba) THEN ! ! Correction Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin ! on modifie le rsol pour que la resistance croisse subitement si on s'approche ! du fond. En gros, rsol=hdry*rsol_cste pour hdry < 1m70 !rsol(ji) = dss(ji,1) * rsol_cste !rsol(ji) = ( drysoil_frac(ji) + un/(10.*(dpu_cste - drysoil_frac(ji))+1.e-10)**2 ) * rsol_cste rsol(ji) = ( hdry(ji) + un/(10.*(dpu_cste - hdry(ji))+1.e-10)**2 ) * rsol_cste ENDIF ENDDO ! ! litter humidity. ! !!$ DO ji = 1, kjpindex !!$ litterhumdiag(ji) = EXP( - dss(ji,1) / hcrit_litter ) !!$ ENDDO litterhumdiag(:) = EXP( - hdry(:) / hcrit_litter ) ! special case: it has just been raining a few drops. The upper soil ! reservoir is ridiculously small, but the dry soil height is zero. ! Don't take it into account. !!$ DO ji = 1, kjpindex !!$ IF ( ( dss(ji,1) .LT. min_sechiba ) .AND. & !!$ ( mean_dsg(ji) .GT. min_sechiba ) .AND. & !!$ ( mean_dsg(ji) .LT. 5.E-4 ) ) THEN !!$ litterhumdiag(ji) = zero !!$ ENDIF !!$ ENDDO WHERE ( ( hdry(:) .LT. min_sechiba ) .AND. & ( mean_dsg(:) .GT. min_sechiba ) .AND. ( mean_dsg(:) .LT. 5.E-4 ) ) litterhumdiag(:) = zero ENDWHERE ! IF (long_print) WRITE (numout,*) ' hydrolc_var_init done ' END SUBROUTINE hydrolc_var_init !! This routine computes snow processes !! SUBROUTINE hydrolc_snow (kjpindex, dtradia, precip_rain, precip_snow , temp_sol_new, soilcap,& & frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, & & tot_melt, snowdepth) ! ! 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) :: precip_rain !! Rainfall REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snow precipitation REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! New soil temperature REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap !! Soil capacity REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ... REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ... ! modified fields REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapnu !! Bare soil evaporation REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow !! Snow mass [Kg/m^2] REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !! Ice water balance REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age!! Snow age on ice, lakes, ... ! output fields REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_melt !! Total melt REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snowdepth !! Snow depth ! ! local declaration ! INTEGER(i_std) :: ji, jv REAL(r_std), DIMENSION (kjpindex) :: d_age !! Snow age change REAL(r_std), DIMENSION (kjpindex) :: xx !! temporary REAL(r_std) :: snowmelt_tmp !! The name says it all ! LOGICAL, DIMENSION (kjpindex) :: warnings LOGICAL :: any_warning ! ! for continental points ! ! ! 0. initialisation ! DO jv = 1, nnobio DO ji=1,kjpindex subsnownobio(ji,jv) = zero ENDDO ENDDO DO ji=1,kjpindex subsnowveg(ji) = zero snowmelt(ji) = zero icemelt(ji) = zero tot_melt(ji) = zero ENDDO ! ! 1. On vegetation ! !cdir NODEP DO ji=1,kjpindex ! ! 1.1. It is snowing ! snow(ji) = snow(ji) + (un - totfrac_nobio(ji))*precip_snow(ji) ENDDO ! DO ji=1,kjpindex ! ! 1.2. Sublimation - separate between vegetated and no-veget fractions ! Care has to be taken as we might have sublimation from the ! the frac_nobio while there is no snow on the rest of the grid. ! IF ( snow(ji) > snowcri ) THEN subsnownobio(ji,iice) = frac_nobio(ji,iice)*vevapsno(ji) subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice) ELSE ! Correction Nathalie - Juillet 2006. ! On doit d'abord tester s'il existe un frac_nobio! ! Pour le moment je ne regarde que le iice IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN subsnownobio(ji,iice) = vevapsno(ji) subsnowveg(ji) = zero ELSE subsnownobio(ji,iice) = zero subsnowveg(ji) = vevapsno(ji) ENDIF ENDIF ENDDO ! warnings(:) = .FALSE. any_warning = .FALSE. !cdir NODEP DO ji=1,kjpindex ! ! 1.2.1 Check that sublimation on the vegetated fraction is possible. ! IF (subsnowveg(ji) .GT. snow(ji)) THEN ! What could not be sublimated goes into bare soil evaporation ! Nathalie - Juillet 2006 - il faut avant tout tester s'il existe du ! frac_nobio sur ce pixel pour eviter de puiser dans le sol! IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN subsnownobio(ji,iice) = subsnownobio(ji,iice) + (subsnowveg(ji) - snow(ji)) ELSE vevapnu(ji) = vevapnu(ji) + (subsnowveg(ji) - snow(ji)) warnings(ji) = .TRUE. any_warning = .TRUE. ENDIF ! Sublimation is thus limited to what is available subsnowveg(ji) = snow(ji) snow(ji) = zero vevapsno(ji) = subsnowveg(ji) + subsnownobio(ji,iice) ELSE snow(ji) = snow(ji) - subsnowveg(ji) ENDIF ENDDO IF ( any_warning ) THEN WRITE(numout,*)' ATTENTION on prend de l eau au sol nu sur au moins un point car evapsno est trop fort!' !!$ DO ji=1,kjpindex !!$ IF ( warnings(ji) ) THEN !!$ WRITE(numout,*)' ATTENTION on prend de l eau au sol nu car evapsno est trop fort!' !!$ WRITE(numout,*)' ',ji,' vevapnu (en mm/jour) = ',vevapnu(ji)*one_day/dtradia !!$ ENDIF !!$ ENDDO ENDIF ! warnings(:) = .FALSE. any_warning = .FALSE. !cdir NODEP DO ji=1,kjpindex ! ! 1.3. snow melt only if temperature positive ! IF (temp_sol_new(ji).GT.tp_00) THEN ! IF (snow(ji).GT.sneige) THEN ! snowmelt(ji) = (un - frac_nobio(ji,iice))*(temp_sol_new(ji) - tp_00) * soilcap(ji) / chalfu0 ! ! 1.3.1.1 enough snow for melting or not ! IF (snowmelt(ji).LT.snow(ji)) THEN snow(ji) = snow(ji) - snowmelt(ji) ELSE snowmelt(ji) = snow(ji) snow(ji) = zero END IF ! ELSEIF (snow(ji).GE.zero) THEN ! ! 1.3.2 not enough snow ! snowmelt(ji) = snow(ji) snow(ji) = zero ELSE ! ! 1.3.3 negative snow - now snow melt ! snow(ji) = zero snowmelt(ji) = zero warnings(ji) = .TRUE. any_warning = .TRUE. ! END IF ENDIF ENDDO IF ( any_warning ) THEN DO ji=1,kjpindex IF ( warnings(ji) ) THEN WRITE(numout,*) 'hydrolc_snow: WARNING! snow was negative and was reset to zero for point ',ji,'. ' ENDIF ENDDO ENDIF ! DO ji=1,kjpindex ! ! 1.4. Ice melt only if there is more than a given mass : maxmass_glacier, ! i.e. only weight melts glaciers ! ! Ajouts Edouard Davin / Nathalie de Noblet add extra to melting ! IF ( snow(ji) .GT. maxmass_glacier ) THEN snowmelt(ji) = snowmelt(ji) + (snow(ji) - maxmass_glacier) snow(ji) = maxmass_glacier ENDIF ! END DO ! ! 2. On Land ice ! DO ji=1,kjpindex ! ! 2.1. It is snowing ! snow_nobio(ji,iice) = snow_nobio(ji,iice) + frac_nobio(ji,iice)*precip_snow(ji) + & & frac_nobio(ji,iice)*precip_rain(ji) ! ! 2.2. Sublimation - was calculated before it can give us negative snow_nobio but that is OK ! Once it goes below a certain values (-maxmass_glacier for instance) we should kill ! the frac_nobio(ji,iice) ! ! snow_nobio(ji,iice) = snow_nobio(ji,iice) - subsnownobio(ji,iice) ! ! 2.3. snow melt only for continental ice fraction ! snowmelt_tmp = zero IF (temp_sol_new(ji) .GT. tp_00) THEN ! ! 2.3.1 If there is snow on the ice-fraction it can melt ! snowmelt_tmp = frac_nobio(ji,iice)*(temp_sol_new(ji) - tp_00) * soilcap(ji) / chalfu0 ! IF ( snowmelt_tmp .GT. snow_nobio(ji,iice) ) THEN snowmelt_tmp = MAX( zero, snow_nobio(ji,iice)) ENDIF snowmelt(ji) = snowmelt(ji) + snowmelt_tmp snow_nobio(ji,iice) = snow_nobio(ji,iice) - snowmelt_tmp ! ENDIF ! ! 2.4 Ice melt only if there is more than a given mass : maxmass_glacier, ! i.e. only weight melts glaciers ! ! IF ( snow_nobio(ji,iice) .GT. maxmass_glacier ) THEN icemelt(ji) = snow_nobio(ji,iice) - maxmass_glacier snow_nobio(ji,iice) = maxmass_glacier ENDIF ! END DO ! ! 3. On other surface types - not done yet ! IF ( nnobio .GT. 1 ) THEN WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW' CALL ipslerr (3,'hydrolc_snow', '', & & 'CANNOT TREAT SNOW ON THESE SURFACE TYPES', '') ENDIF ! ! 4. computes total melt (snow and ice) ! DO ji = 1, kjpindex tot_melt(ji) = icemelt(ji) + snowmelt(ji) ENDDO ! ! 5. computes snow age on veg and ice (for albedo) ! DO ji = 1, kjpindex ! ! 5.1 Snow age on vegetation ! IF (snow(ji) .LE. zero) THEN snow_age(ji) = zero ELSE snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dtradia/one_day) & & * EXP(-precip_snow(ji) / snow_trans) ENDIF ! ! 5.2 Snow age on ice ! ! age of snow on ice: a little bit different because in cold regions, we really ! cannot negect the effect of cold temperatures on snow metamorphism any more. ! IF (snow_nobio(ji,iice) .LE. zero) THEN snow_nobio_age(ji,iice) = zero ELSE ! d_age(ji) = ( snow_nobio_age(ji,iice) + & & (un - snow_nobio_age(ji,iice)/max_snow_age) * dtradia/one_day ) * & & EXP(-precip_snow(ji) / snow_trans) - snow_nobio_age(ji,iice) IF (d_age(ji) .GT. zero ) THEN xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero ) xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std d_age(ji) = d_age(ji) / (un+xx(ji)) ENDIF snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero ) ! ENDIF ENDDO ! ! 6.0 Diagnose the depth of the snow layer ! DO ji = 1, kjpindex snowdepth(ji) = snow(ji) /sn_dens ENDDO IF (long_print) WRITE (numout,*) ' hydrolc_snow done ' END SUBROUTINE hydrolc_snow !! This routine computes canopy processes !! SUBROUTINE hydrolc_canop (kjpindex, precip_rain, vevapwet, veget, qsintmax, qsintveg, precisol) ! ! interface description ! INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rain precipitation REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: vevapwet !! Interception loss REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation type REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: qsintmax !! Maximum water on vegetation for interception ! modified fields REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: qsintveg !! Water on vegetation due to interception ! output fields REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: precisol !! Eau tombee sur le sol ! ! local declaration ! INTEGER(i_std) :: ji, jv REAL(r_std), DIMENSION (kjpindex,nvm) :: zqsintvegnew ! calcul de qsintmax a prevoir a chaque pas de temps ! dans ini_sechiba ! boucle sur les points continentaux ! calcul de qsintveg au pas de temps suivant ! par ajout du flux interception loss ! calcule par enerbil en fonction ! des calculs faits dans diffuco ! calcul de ce qui tombe sur le sol ! avec accumulation dans precisol ! essayer d'harmoniser le traitement du sol nu ! avec celui des differents types de vegetation ! fait si on impose qsintmax ( ,1) = 0.0 ! ! loop for continental subdomain ! ! ! 1. evaporation off the continents ! ! 1.1 The interception loss is take off the canopy. DO jv=1,nvm qsintveg(:,jv) = qsintveg(:,jv) - vevapwet(:,jv) END DO ! 1.2 It is raining : precip_rain is shared for each vegetation ! type ! sum (veget (1,nvm)) must be egal to 1-totfrac_nobio. ! iniveget computes veget each day ! DO jv=1,nvm ! Correction Nathalie - Juin 2006 - une partie de la pluie arrivera toujours sur le sol ! sorte de throughfall supplementaire !qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * precip_rain(:) qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * ((1-throughfall_by_pft(jv))*precip_rain(:)) END DO ! ! 1.3 Limits the effect and sum what receives soil ! precisol(:,:) = zero DO jv=1,nvm DO ji = 1, kjpindex zqsintvegnew(ji,jv) = MIN (qsintveg(ji,jv),qsintmax(ji,jv)) ! correction throughfall Nathalie - Juin 2006 !precisol(ji,jv) = qsintveg(ji,jv ) - zqsintvegnew (ji,jv) precisol(ji,jv) = (veget(ji,jv)*throughfall_by_pft(jv)*precip_rain(ji)) + qsintveg(ji,jv ) - zqsintvegnew (ji,jv) ENDDO ENDDO ! ! 1.4 swap qsintveg to the new value ! DO jv=1,nvm qsintveg(:,jv) = zqsintvegnew (:,jv) END DO IF (long_print) WRITE (numout,*) ' hydrolc_canop done ' END SUBROUTINE hydrolc_canop !! !! !! SUBROUTINE hydrolc_vegupd(kjpindex, veget, ruu_ch, qsintveg, gqsb, bqsb, dsg, dss, dsp, resdist) ! ! ! The vegetation cover has changed and we need to adapt the reservoir distribution. ! You may note that this occurs after evaporation and so on have been computed. It is ! not a problem as a new vegetation fraction will start with humrel=0 and thus will have no ! evaporation. If this is not the case it should have been caught above. ! ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex ! input fields REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in) :: veget !! New vegetation map REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: ruu_ch !! Quantite d'eau maximum ! modified fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintveg !! Water on vegetation due to interception REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: gqsb !! Hauteur d'eau dans le reservoir de surface REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: bqsb !! Hauteur d'eau dans le reservoir profond REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsg !! Hauteur du reservoir de surface REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dss !! Hauteur au dessus du reservoir de surface REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsp !! Hauteur au dessus du reservoir profond REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(inout) :: resdist !! Old vegetation map ! ! ! local declaration ! INTEGER(i_std) :: ji,jv ! REAL(r_std),DIMENSION (kjpindex,nvm) :: qsintveg2 !! Water on vegetation due to interception over old veget REAL(r_std), DIMENSION (kjpindex,nvm) :: bdq, gdq, qsdq REAL(r_std), DIMENSION (kjpindex,nvm) :: vmr !! variation of veget REAL(r_std), DIMENSION(kjpindex) :: gtr, btr, vtr, qstr, fra REAL(r_std), DIMENSION(kjpindex) :: vegchtot REAL(r_std), PARAMETER :: EPS1 = EPSILON(un) ! ! DO jv = 1, nvm DO ji = 1, kjpindex IF ( ABS(veget(ji,jv)-resdist(ji,jv)) .GT. EPS1 ) THEN vmr(ji,jv) = veget(ji,jv)-resdist(ji,jv) ELSE vmr(ji,jv) = zero ENDIF ! IF (resdist(ji,jv) .GT. zero) THEN qsintveg2(ji,jv) = qsintveg(ji,jv)/resdist(ji,jv) ELSE qsintveg2(ji,jv) = zero ENDIF ENDDO ENDDO ! vegchtot(:) = zero DO jv = 1, nvm DO ji = 1, kjpindex vegchtot(ji) = vegchtot(ji) + ABS( vmr(ji,jv) ) ENDDO ENDDO ! DO jv = 1, nvm DO ji = 1, kjpindex IF ( vegchtot(ji) .GT. zero ) THEN gdq(ji,jv) = ABS(vmr(ji,jv)) * gqsb(ji,jv) bdq(ji,jv) = ABS(vmr(ji,jv)) * bqsb(ji,jv) qsdq(ji,jv) = ABS(vmr(ji,jv)) * qsintveg2(ji,jv) ENDIF ENDDO ENDDO ! ! calculate water mass that we have to redistribute ! gtr(:) = zero btr(:) = zero qstr(:) = zero vtr(:) = zero ! ! DO jv = 1, nvm DO ji = 1, kjpindex IF ( ( vegchtot(ji) .GT. zero ) .AND. ( vmr(ji,jv) .LT. zero ) ) THEN gtr(ji) = gtr(ji) + gdq(ji,jv) btr(ji) = btr(ji) + bdq(ji,jv) qstr(ji) = qstr(ji) + qsdq(ji,jv) vtr(ji) = vtr(ji) - vmr(ji,jv) ENDIF ENDDO ENDDO ! ! put it into reservoir of plant whose surface area has grown DO jv = 1, nvm DO ji = 1, kjpindex IF ( vegchtot(ji) .GT. zero .AND. ABS(vtr(ji)) .GT. EPS1) THEN fra(ji) = vmr(ji,jv) / vtr(ji) IF ( vmr(ji,jv) .GT. zero) THEN IF (veget(ji,jv) .GT. zero) THEN gqsb(ji,jv) = (resdist(ji,jv)*gqsb(ji,jv) + fra(ji)*gtr(ji))/veget(ji,jv) bqsb(ji,jv) = (resdist(ji,jv)*bqsb(ji,jv) + fra(ji)*btr(ji))/veget(ji,jv) ENDIF qsintveg(ji,jv) = qsintveg(ji,jv) + fra(ji)* qstr(ji) ELSE qsintveg(ji,jv) = qsintveg(ji,jv) - qsdq(ji,jv) ENDIF ! ! dss is not altered so that this transfer of moisture does not directly ! affect transpiration IF (gqsb(ji,jv) .LT. min_sechiba) THEN dsg(ji,jv) = zero ELSE dsg(ji,jv) = (dss(ji,jv) * ruu_ch(ji) + gqsb(ji,jv)) & / ruu_ch(ji) ENDIF dsp(ji,jv) = dpu_cste - bqsb(ji,jv) / ruu_ch(ji) ENDIF ENDDO ENDDO ! Now that the work is done resdist needs an update ! DO jv = 1, nvm resdist(:,jv) = veget(:,jv) ENDDO ! ! Where vegetation fraction is zero, set water to that of bare soil. ! This does not create any additional water. ! DO jv = 2, nvm DO ji = 1, kjpindex IF ( veget(ji,jv) .LT. EPS1 ) THEN gqsb(ji,jv) = gqsb(ji,1) bqsb(ji,jv) = bqsb(ji,1) dsg(ji,jv) = dsg(ji,1) dss(ji,jv) = dss(ji,1) dsp(ji,jv) = dsp(ji,1) ENDIF ENDDO ENDDO RETURN ! END SUBROUTINE hydrolc_vegupd !! !! This routines computes soil processes !! SUBROUTINE hydrolc_soil(kjpindex, vevapnu, precisol, returnflow, irrigation, tot_melt, mx_eau_var, veget, ruu_ch, transpir,& & gqsb, bqsb, dsg, dss, rsol, drysoil_frac, hdry, dsp, runoff, run_off_tot, drainage, humrel, vegstress, & & shumdiag, litterhumdiag) ! ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex ! input fields REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: vevapnu !! Bare soil evaporation REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: precisol !! Eau tombee sur le sol REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: returnflow !! Water returning to the deep reservoir REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: irrigation !! Irrigation REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_melt !! Total melt REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: mx_eau_var !! Profondeur du reservoir contenant le !! maximum d'eau REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in) :: veget !! Vegetation map REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: ruu_ch !! Quantite d'eau maximum REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: transpir !! Transpiration ! modified fields REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: gqsb !! Hauteur d'eau dans le reservoir de surface REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: bqsb !! Hauteur d'eau dans le reservoir profond REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsg !! Hauteur du reservoir de surface REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dss !! Hauteur au dessus du reservoir de surface REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsp !! Hauteur au dessus du reservoir profond ! output fields REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: runoff !! Ruissellement REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: run_off_tot !! Complete runoff REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: drainage !! Drainage REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out) :: humrel !! Relative humidity REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out) :: vegstress !! Veg. moisture stress (only for vegetation growth) REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: shumdiag !! relative soil moisture REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: rsol !! resistance to bare soil evaporation REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: drysoil_frac !! Fraction of visible fry soil REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: hdry !! Dry soil heigth in meters ! ! local declaration ! INTEGER(i_std) :: ji,jv, jd REAL(r_std), DIMENSION(kjpindex) :: zhumrel_lo, zhumrel_up REAL(r_std), DIMENSION(kjpindex,nvm) :: zeflux !! Soil evaporation REAL(r_std), DIMENSION(kjpindex,nvm) :: zpreci !! Soil precipitation LOGICAL, DIMENSION(kjpindex,nvm) :: warning REAL(r_std) :: gtr, btr REAL(r_std), DIMENSION(kjpindex) :: mean_dsg LOGICAL :: OnceMore INTEGER(i_std), PARAMETER :: nitermax = 100 INTEGER(i_std) :: niter INTEGER(i_std) :: nbad LOGICAL, DIMENSION(kjpindex,nvm) :: lbad LOGICAL, DIMENSION(kjpindex) :: lbad_ij REAL(r_std) :: gqseuil , eausup, wd1 REAL(r_std), DIMENSION(nbdl+1) :: tmp_dl ! Ajout Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin ! Modifs stabilite REAL(r_std), DIMENSION(kjpindex,nvm) :: a_subgrd ! ! 0. we have only one flux field corresponding to water evaporated from the surface ! DO jv=1,nvm DO ji = 1, kjpindex IF ( veget(ji,jv) .GT. zero ) THEN zeflux(ji,jv) = transpir(ji,jv)/veget(ji,jv) zpreci(ji,jv) = precisol(ji,jv)/veget(ji,jv) ELSE zeflux(ji,jv) = zero zpreci(ji,jv) = zero ENDIF ENDDO ENDDO ! ! We need a test on the bare soil fraction because we can have bare soil evaporation even when ! there is no bare soil because of transfers (snow for instance). This should only apply if there ! is vegetation but we do not test this case. ! ! case 1 - there is vegetation and bare soil DO ji = 1, kjpindex IF ( (vegtot(ji) .GT. zero) .AND. (veget(ji,1) .GT. min_sechiba) ) THEN zeflux(ji,1) = vevapnu(ji)/veget(ji,1) ENDIF ENDDO ! case 2 - there is vegetation but no bare soil DO jv = 1, nvm DO ji = 1, kjpindex IF ( (vegtot(ji) .GT. zero) .AND. (veget(ji,1) .LE. min_sechiba)& & .AND. (veget(ji,jv) .GT. min_sechiba)) THEN zeflux(ji,jv) = zeflux(ji,jv) + vevapnu(ji)/vegtot(ji) ENDIF ENDDO ENDDO ! ! 0.1 Other temporary variables ! tmp_dl(1) = 0 tmp_dl(2:nbdl+1) = diaglev(1:nbdl) ! ! ! 1.1 Transpiration for each vegetation type ! Evaporated water is taken out of the ground. ! ! DO jv=1,nvm DO ji=1,kjpindex ! gqsb(ji,jv) = gqsb(ji,jv) - zeflux(ji,jv) ! ! 1.2 Add snow and ice melt, troughfall from vegetation and irrigation. ! IF(vegtot(ji) .NE. zero) THEN ! snow and ice melt and troughfall from vegetation gqsb(ji,jv) = gqsb(ji,jv) + zpreci(ji,jv) + tot_melt(ji)/vegtot(ji) ! ! We take care to add the irrigation only to the vegetated part if possible ! IF (ABS(vegtot(ji)-veget(ji,1)) .LE. min_sechiba) THEN gqsb(ji,jv) = gqsb(ji,jv) + irrigation(ji)/vegtot(ji) ELSE IF ( jv > 1 ) THEN ! Only add the irrigation to the upper soil if there is a reservoir. ! Without this the water evaporates right away. IF ( gqsb(ji,jv) > zero ) THEN gqsb(ji,jv) = gqsb(ji,jv) + irrigation(ji)/(vegtot(ji)-veget(ji,1)) ELSE bqsb(ji,jv) = bqsb(ji,jv) + irrigation(ji)/(vegtot(ji)-veget(ji,1)) ENDIF ENDIF ENDIF ! ! 1.3 We add the water returned from rivers to the lower reservoir. ! bqsb(ji,jv) = bqsb(ji,jv) + returnflow(ji)/vegtot(ji) ! ENDIF ! END DO ENDDO ! ! 1.3 Computes run-off ! runoff(:,:) = zero ! ! 1.4 Soil moisture is updated ! warning(:,:) = .FALSE. DO jv=1,nvm DO ji = 1, kjpindex ! runoff(ji,jv) = MAX(gqsb(ji,jv) + bqsb(ji,jv) - mx_eau_var(ji), zero) ! IF (mx_eau_var(ji) .LE. (gqsb(ji,jv) + bqsb(ji,jv))) THEN ! ! 1.4.1 Plus de reservoir de surface: le sol est sature ! d'eau. Le reservoir de surface est inexistant ! Tout est dans le reservoir de fond. ! Le ruissellement correspond a ce qui deborde. ! gqsb(ji,jv) = zero dsg(ji,jv) = zero bqsb(ji,jv) = mx_eau_var(ji) dsp(ji,jv) = zero dss(ji,jv) = dsp (ji,jv) ELSEIF ((gqsb(ji,jv) + bqsb(ji,jv)).GE.zero) THEN ! IF (gqsb(ji,jv) .GT. dsg(ji,jv) * ruu_ch(ji)) THEN ! ! 1.4.2 On agrandit le reservoir de surface ! car il n y a pas eu ruissellement ! et toute l'eau du reservoir de surface ! tient dans un reservoir dont la taille ! est plus importante que l actuel. ! La hauteur de sol sec dans le reservoir ! de surface est alors nulle. ! Le reste ne bouge pas. ! dsg(ji,jv) = gqsb(ji,jv) / ruu_ch(ji) dss(ji,jv) = zero ELSEIF (gqsb(ji,jv) .GT. zero ) THEN ! ! 1.4.3 L eau tient dans le reservoir de surface ! tel qu il existe. ! Calcul de la nouvelle hauteur de sol sec ! dans la couche de surface. ! Le reste ne bouge pas. ! dss(ji,jv) = ((dsg(ji,jv) * ruu_ch(ji)) - gqsb(ji,jv)) / ruu_ch(ji) ELSE ! ! 1.4.4 La quantite d eau du reservoir de surface ! est negative. Cela revient a enlever ! cette quantite au reservoir profond. ! Le reservoir de surface est alors vide. ! (voir aussi 1.4.1) ! bqsb(ji,jv) = bqsb(ji,jv) + gqsb(ji,jv) dsp(ji,jv) = dpu_cste - bqsb(ji,jv) / ruu_ch(ji) gqsb(ji,jv) = zero dsg(ji,jv) = zero dss(ji,jv) = dsp(ji,jv) END IF ELSE ! ! 1.4.5 Le reservoir profond est aussi asseche. ! La quantite d eau a enlever depasse la quantite ! disponible dans le reservoir profond. ! ! ! Ceci ne devrait jamais arriver plus d'une fois par point. C-a-d une fois la valeur negative ! atteinte les flux doivent etre nuls. On ne signal que ce cas donc. ! IF ( ( zeflux(ji,jv) .GT. zero ) .AND. & ( gqsb(ji,jv) + bqsb(ji,jv) .LT. -1.e-10 ) ) THEN warning(ji,jv) = .TRUE. ! WRITE (numout,*) 'WARNING! Soil Moisture will be negative' ! WRITE (numout,*) 'ji, jv = ', ji,jv ! WRITE (numout,*) 'mx_eau_var = ', mx_eau_var(ji) ! WRITE (numout,*) 'veget, resdist =', veget(ji,jv), resdist(ji,jv) ! WRITE (numout,*) 'bqsb = ', bqsb(ji,jv) ! WRITE (numout,*) 'gqsb = ', gqsb(ji,jv) ! WRITE (numout,*) 'dss = ', dss(ji,jv) ! WRITE (numout,*) 'dsg = ', dsg(ji,jv) ! WRITE (numout,*) 'dsp = ', dsp(ji,jv) ! WRITE (numout,*) 'humrel = ', humrel(ji, jv) ! WRITE (numout,*) 'Soil evaporation = ', zeflux(ji,jv) ! WRITE (numout,*) 'input = ',precisol(ji, jv), tot_melt(ji) ! WRITE (numout,*) '============================' ENDIF ! bqsb(ji,jv) = gqsb(ji,jv) + bqsb(ji,jv) dsp(ji,jv) = dpu_cste gqsb(ji,jv) = zero dsg(ji,jv) = zero dss(ji,jv) = dsp(ji,jv) ! ENDIF ! ENDDO ENDDO ! nbad = COUNT( warning(:,:) .EQV. .TRUE. ) IF ( nbad .GT. 0 ) THEN WRITE(numout,*) 'hydrolc_soil: WARNING! Soil moisture was negative at', & nbad, ' points:' !DO jv = 1, nvm ! DO ji = 1, kjpindex ! IF ( warning(ji,jv) ) THEN ! WRITE(numout,*) ' ji,jv = ', ji,jv ! WRITE (numout,*) 'mx_eau_var = ', mx_eau_var(ji) ! WRITE (numout,*) 'veget, resdist =', veget(ji,jv), resdist(ji,jv) ! WRITE (numout,*) 'bqsb = ', bqsb(ji,jv) ! WRITE (numout,*) 'gqsb = ', gqsb(ji,jv) ! WRITE (numout,*) 'runoff = ',runoff(ji,jv) ! WRITE (numout,*) 'dss = ', dss(ji,jv) ! WRITE (numout,*) 'dsg = ', dsg(ji,jv) ! WRITE (numout,*) 'dsp = ', dsp(ji,jv) ! WRITE (numout,*) 'humrel = ', humrel(ji, jv) ! WRITE (numout,*) 'Soil evaporation = ', zeflux(ji,jv) ! WRITE (numout,*) 'Soil precipitation = ',zpreci(ji,jv) ! WRITE (numout,*) 'input = ',precisol(ji, jv), tot_melt(ji) ! WRITE (numout,*) 'returnflow = ',returnflow(ji) ! WRITE (numout,*) '============================' ! ENDIF ! ENDDO !ENDDO ENDIF ! ! 2.0 Very large upper reservoirs or very close upper and lower reservoirs ! can be deadlock situations for Choisnel. They are handled here ! IF (long_print) WRITE(numout,*) 'hydro_soil 2.0 : Resolve deadlocks' DO jv=1,nvm DO ji=1,kjpindex ! ! 2.1 the two reservoirs are very close to each other ! IF ( ABS(dsp(ji,jv)-dsg(ji,jv)) .LT. min_resdis ) THEN bqsb(ji,jv) = bqsb(ji,jv) + gqsb(ji,jv) dsp(ji,jv) = dpu_cste - bqsb(ji,jv) / ruu_ch(ji) gqsb(ji,jv) = zero dsg(ji,jv) = zero dss(ji,jv) = dsp(ji,jv) ENDIF ! ! 2.2 Draine some water from the upper to the lower reservoir ! gqseuil = min_resdis * deux * ruu_ch(ji) eausup = dsg(ji,jv) * ruu_ch(ji) wd1 = .75*eausup ! IF (eausup .GT. gqseuil) THEN gdrainage(ji,jv) = min_drain * (gqsb(ji,jv)/eausup) ! IF ( gqsb(ji,jv) .GE. wd1 .AND. dsg(ji,jv) .GT. 0.10 ) THEN ! gdrainage(ji,jv) = gdrainage(ji,jv) + & (max_drain-min_drain)*((gqsb(ji,jv)-wd1) / (eausup-wd1))**exp_drain ! ENDIF ! gdrainage(ji,jv)=MIN(gdrainage(ji,jv), MAX(gqsb(ji,jv), zero)) ! ELSE gdrainage(ji,jv)=zero ENDIF ! gqsb(ji,jv) = gqsb(ji,jv) - gdrainage(ji,jv) bqsb(ji,jv) = bqsb(ji,jv) + gdrainage(ji,jv) dsg(ji,jv) = dsg(ji,jv) - gdrainage(ji,jv) / ruu_ch(ji) dsp(ji,jv) = dpu_cste - bqsb(ji,jv)/ruu_ch(ji) ! ! ENDDO ! ENDDO ! ! 3.0 Diffusion of water between the reservoirs of the different plants ! IF (long_print) WRITE(numout,*) 'hydrolc_soil 3.0 : Vertical diffusion' mean_bqsb(:) = zero mean_gqsb(:) = zero DO jv = 1, nvm DO ji = 1, kjpindex IF ( vegtot(ji) .GT. zero ) THEN mean_bqsb(ji) = mean_bqsb(ji) + veget(ji,jv)/vegtot(ji)*bqsb(ji,jv) mean_gqsb(ji) = mean_gqsb(ji) + veget(ji,jv)/vegtot(ji)*gqsb(ji,jv) ENDIF ENDDO ENDDO OnceMore = .TRUE. niter = 0 !YM lbad_ij(:)=.TRUE. ! nitermax prevents infinite loops (should actually never occur) DO WHILE ( OnceMore .AND. ( niter .LT. nitermax ) ) ! niter = niter + 1 !YM DO jv = 1, nvm ! DO ji = 1, kjpindex IF (lbad_ij(ji)) THEN IF ( veget(ji,jv) .GT. zero ) THEN ! bqsb(ji,jv) = mean_bqsb(ji) dsp(ji,jv) = dpu_cste - bqsb(ji,jv)/ruu_ch(ji) ENDIF ENDIF ! ENDDO ENDDO ! ! where do we have to do something? lbad(:,:) = ( ( dsp(:,:) .LT. dsg(:,:) ) .AND. & ( dsg(:,:) .GT. zero ) .AND. & ( veget(:,:) .GT. zero ) ) ! ! if there are no such points any more, we'll do no more iteration IF ( COUNT( lbad(:,:) ) .EQ. 0 ) OnceMore = .FALSE. DO ji = 1, kjpindex IF (COUNT(lbad(ji,:)) == 0 ) lbad_ij(ji)=.FALSE. ENDDO ! DO jv = 1, nvm !YM ! ! ! DO ji = 1, kjpindex ! IF ( veget(ji,jv) .GT. zero ) THEN ! ! ! bqsb(ji,jv) = mean_bqsb(ji) ! dsp(ji,jv) = dpu_cste - bqsb(ji,jv)/ruu_ch(ji) ! ENDIF ! ! ! ENDDO ! DO ji = 1, kjpindex IF ( lbad(ji,jv) ) THEN ! runoff(ji,jv) = runoff(ji,jv) + & MAX( bqsb(ji,jv) + gqsb(ji,jv) - mx_eau_var(ji), zero) ! bqsb(ji,jv) = MIN( bqsb(ji,jv) + gqsb(ji,jv), mx_eau_var(ji)) ! gqsb(ji,jv) = zero dsp(ji,jv) = dpu_cste - bqsb(ji,jv)/ruu_ch(ji) dss(ji,jv) = dsp(ji,jv) dsg(ji,jv) = zero ! ENDIF ENDDO ! ENDDO ! mean_bqsb(:) = zero mean_gqsb(:) = zero DO jv = 1, nvm DO ji = 1, kjpindex IF ( vegtot(ji) .GT. zero ) THEN mean_bqsb(ji) = mean_bqsb(ji) + veget(ji,jv)/vegtot(ji)*bqsb(ji,jv) mean_gqsb(ji) = mean_gqsb(ji) + veget(ji,jv)/vegtot(ji)*gqsb(ji,jv) ENDIF ENDDO ENDDO ! ENDDO ! ! 4. computes total runoff ! IF (long_print) WRITE(numout,*) 'hydrolc_soil 4.0: Computes total runoff' run_off_tot(:) = zero DO ji = 1, kjpindex IF ( vegtot(ji) .GT. zero ) THEN DO jv = 1, nvm run_off_tot(ji) = run_off_tot(ji) + (runoff(ji,jv)*veget(ji,jv)) ENDDO ELSE run_off_tot(ji) = tot_melt(ji) + irrigation(ji) ENDIF ENDDO ! ! 4.1 We estimate some drainage ! ! drainage(:) = 0.95 * run_off_tot(:) run_off_tot(:) = run_off_tot(:) - drainage(:) ! ! 5. Some averaged diagnostics ! IF (long_print) WRITE(numout,*) 'hydro_soil 5.0: Diagnostics' ! ! 5.1 reset dsg if necessary ! WHERE (gqsb(:,:) .LE. zero) dsg(:,:) = zero ! DO ji=1,kjpindex ! ! 5.2 Compute an average moisture profile ! mean_dsg(ji) = mean_gqsb(ji)/ruu_ch(ji) ! ENDDO ! ! 6. Compute the moisture stress on vegetation ! IF (long_print) WRITE(numout,*) 'hydro_soil 6.0 : Moisture stress' a_subgrd(:,:) = zero DO jv = 1, nvm DO ji=1,kjpindex ! ! computes relative surface humidity ! ! Only use the standard formulas if total soil moisture is larger than zero. ! Else stress functions are set to zero. ! This will avoid that large negative soil moisture accumulate over time by the ! the creation of small skin reservoirs which evaporate quickly. ! IF ( gqsb(ji,jv)+bqsb(ji,jv) .GT. zero ) THEN ! IF (dsg(ji,jv).EQ. zero .OR. gqsb(ji,jv).EQ.zero) THEN humrel(ji,jv) = EXP( - humcste(jv) * dpu_cste * (dsp(ji,jv)/dpu_cste) ) dsg(ji,jv) = zero ! ! if the dry soil height is larger than the one corresponding ! to the wilting point, or negative lower soil moisture : humrel is 0.0 ! IF (dsp(ji,jv).GT.(dpu_cste - (qwilt / ruu_ch(ji))) .OR. bqsb(ji,jv).LT.zero) THEN humrel(ji,jv) = zero ENDIF ! ! In this case we can take for vegetation growth the same values for humrel and vegstress ! vegstress(ji,jv) = humrel(ji,jv) ! ELSE ! Corrections Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin !zhumrel_lo(ji) = EXP( - humcste(jv) * dpu_cste * (dsp(ji,jv)/dpu_cste) ) !zhumrel_up(ji) = EXP( - humcste(jv) * dpu_cste * (dss(ji,jv)/dsg(ji,jv)) ) !humrel(ji,jv) = MAX(zhumrel_lo(ji),zhumrel_up(ji)) ! ! As we need a slower variable for vegetation growth the stress is computed ! differently than in humrel. ! zhumrel_lo(ji) = EXP( - humcste(jv) * dsp(ji,jv)) zhumrel_up(ji) = EXP( - humcste(jv) * dss(ji,jv)) ! Ajouts Nathalie - Fred - le 28 Mars 2006 a_subgrd(ji,jv)=MIN(MAX(dsg(ji,jv)-dss(ji,jv),zero)/dsg_min,un) humrel(ji,jv)=a_subgrd(ji,jv)*zhumrel_up(ji)+(un-a_subgrd(ji,jv))*zhumrel_lo(ji) ! vegstress(ji,jv) = zhumrel_lo(ji) + zhumrel_up(ji) - EXP( - humcste(jv) * dsg(ji,jv) ) ! ENDIF ! ELSE ! humrel(ji,jv) = zero vegstress(ji,jv) = zero ! ENDIF ! ENDDO ENDDO ! ! 7. Diagnostics which are needed to carry information to other modules ! ! 7.2 Relative soil moisture ! DO jd = 1,nbdl DO ji = 1, kjpindex IF ( tmp_dl(jd+1) .LT. mean_dsg(ji)) THEN shumdiag(ji,jd) = mean_gqsb(ji)/mx_eau_var(ji) ELSE IF ( tmp_dl(jd) .LT. mean_dsg(ji)) THEN gtr = (mean_dsg(ji)-tmp_dl(jd))/(tmp_dl(jd+1)-tmp_dl(jd)) btr = 1 - gtr shumdiag(ji,jd) = gtr*mean_gqsb(ji)/mx_eau_var(ji) + & & btr*mean_bqsb(ji)/mx_eau_var(ji) ELSE shumdiag(ji,jd) = mean_bqsb(ji)/mx_eau_var(ji) ENDIF ENDIF shumdiag(ji,jd) = MAX(MIN(shumdiag(ji,jd), un), zero) ENDDO ENDDO ! The fraction of visibly dry soil (dry when dss(:,1) = 0.1 m) drysoil_frac(:) = MIN(MAX(dss(:,1),zero)*10._r_std, un) ! Correction Nathalie - le 28 Mars 2006 - re-ecriture drysoil_frac/hdry - Fred Hourdin ! revu 22 novembre 2007 hdry(:) = a_subgrd(:,1)*dss(:,1) + (un-a_subgrd(:,1))*dsp(:,1) ! ! Compute the resistance to bare soil evaporation. ! rsol(:) = -un DO ji = 1, kjpindex IF (veget(ji,1) .GE. min_sechiba) THEN ! ! Correction Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin ! on modifie le rsol pour que la resistance croisse subitement si on s'approche ! du fond. En gros, rsol=hdry*rsol_cste pour hdry < 1m70 !rsol(ji) = dss(ji,1) * rsol_cste rsol(ji) = ( hdry(ji) + un/(10.*(dpu_cste - hdry(ji))+1.e-10)**2 ) * rsol_cste ENDIF ENDDO ! ! 8. litter humidity. ! litterhumdiag(:) = EXP( - hdry(:) / hcrit_litter ) ! special case: it has just been raining a few drops. The upper soil ! reservoir is ridiculously small, but the dry soil height is zero. ! Don't take it into account. WHERE ( ( hdry(:) .LT. min_sechiba ) .AND. & ( mean_dsg(:) .GT. min_sechiba ) .AND. ( mean_dsg(:) .LT. 5.E-4 ) ) litterhumdiag(:) = zero ENDWHERE ! IF (long_print) WRITE (numout,*) ' hydrolc_soil done ' END SUBROUTINE hydrolc_soil !! !! This routines checks the water balance. First it gets the total !! amount of water and then it compares the increments with the fluxes. !! The computation is only done over the soil area as over glaciers (and lakes?) !! we do not have water conservation. !! !! This verification does not make much sense in REAL*4 as the precision is the same as some !! of the fluxes !! SUBROUTINE hydrolc_waterbal (kjpindex, index, first_call, dtradia, veget, totfrac_nobio, qsintveg, snow, snow_nobio,& & precip_rain, precip_snow, returnflow, irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno, run_off_tot, & & drainage) ! ! ! INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map LOGICAL, INTENT (in) :: first_call !! At which time is this routine called ? REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds ! REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio!! Total fraction of continental ice+lakes+... REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass [Kg/m^2] REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !!Ice water balance ! REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: returnflow !! Water returning from routing to the deep reservoir REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: irrigation !! Water from irrigation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tot_melt !! Total melt ! REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vevapwet !! Interception loss REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapnu !! Bare soil evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapsno !! Snow evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: run_off_tot !! Total runoff REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: drainage !! Drainage ! ! LOCAL ! INTEGER(i_std) :: ji, jv, jn REAL(r_std) :: allowed_err REAL(r_std),DIMENSION (kjpindex) :: watveg, delta_water, tot_flux, sum_snow_nobio, sum_vevapwet, sum_transpir ! ! ! IF ( first_call ) THEN tot_water_beg(:) = zero watveg(:) = zero sum_snow_nobio(:) = zero !cdir NODEP DO jv = 1, nvm watveg(:) = watveg(:) + qsintveg(:,jv) ENDDO !cdir NODEP DO jn = 1, nnobio sum_snow_nobio(:) = sum_snow_nobio(:) + snow_nobio(:,jn) ENDDO !cdir NODEP DO ji = 1, kjpindex tot_water_beg(ji) = (mean_bqsb(ji) + mean_gqsb(ji))*vegtot(ji) + & & watveg(ji) + snow(ji) + sum_snow_nobio(ji) ENDDO tot_water_end(:) = tot_water_beg(:) RETURN ENDIF ! ! Check the water balance ! tot_water_end(:) = zero ! DO ji = 1, kjpindex ! ! If the fraction of ice, lakes, etc. does not complement the vegetation fraction then we do not ! need to go any further ! ! Modif Nathalie ! IF ( (un - (totfrac_nobio(ji) + vegtot(ji))) .GT. EPSILON(un) ) THEN IF ( (un - (totfrac_nobio(ji) + vegtot(ji))) .GT. (100*EPSILON(un)) ) THEN WRITE(numout,*) 'HYDROL problem in vegetation or frac_nobio on point ', ji WRITE(numout,*) 'totfrac_nobio : ', totfrac_nobio(ji) WRITE(numout,*) 'vegetation fraction : ', vegtot(ji) !STOP 'in hydrolc_waterbal' ENDIF ENDDO ! watveg(:) = zero sum_vevapwet(:) = zero sum_transpir(:) = zero sum_snow_nobio(:) = zero !cdir NODEP DO jv = 1,nvm watveg(:) = watveg(:) + qsintveg(:,jv) sum_vevapwet(:) = sum_vevapwet(:) + vevapwet(:,jv) sum_transpir(:) = sum_transpir(:) + transpir(:,jv) ENDDO !cdir NODEP DO jn = 1,nnobio sum_snow_nobio(:) = sum_snow_nobio(:) + snow_nobio(:,jn) ENDDO ! !cdir NODEP DO ji = 1, kjpindex tot_water_end(ji) = (mean_bqsb(ji) + mean_gqsb(ji))*vegtot(ji) + & & watveg(ji) + snow(ji) + sum_snow_nobio(ji) ENDDO ! DO ji = 1, kjpindex ! delta_water(ji) = tot_water_end(ji) - tot_water_beg(ji) ! tot_flux(ji) = precip_rain(ji) + precip_snow(ji) + returnflow(ji) + irrigation(ji) - & & sum_vevapwet(ji) - sum_transpir(ji) - vevapnu(ji) - vevapsno(ji) - & & run_off_tot(ji) - drainage(ji) ! ENDDO ! ! Set some precision ! This is a wild guess and corresponds to what works on an IEEE machine ! under double precision (REAL*8). ! allowed_err = 50000*EPSILON(un) ! DO ji = 1, kjpindex IF ( ABS(delta_water(ji)-tot_flux(ji)) .GT. allowed_err ) THEN WRITE(numout,*) 'HYDROL does not conserve water. The erroneous point is : ', ji WRITE(numout,*) 'The error in mm/d is :', (delta_water(ji)-tot_flux(ji))/dtradia*one_day, & & ' and in mm/dt : ', delta_water(ji)-tot_flux(ji) WRITE(numout,*) 'delta_water : ', delta_water(ji), ' tot_flux : ', tot_flux(ji) WRITE(numout,*) 'Actual and allowed error : ', ABS(delta_water(ji)-tot_flux(ji)), allowed_err WRITE(numout,*) 'vegtot : ', vegtot(ji) WRITE(numout,*) 'precip_rain : ', precip_rain(ji) WRITE(numout,*) 'precip_snow : ', precip_snow(ji) WRITE(numout,*) 'Water from irrigation : ', returnflow(ji), irrigation(ji) WRITE(numout,*) 'Total water in soil :', mean_bqsb(ji) + mean_gqsb(ji) WRITE(numout,*) 'Water on vegetation :', watveg(ji) WRITE(numout,*) 'Snow mass :', snow(ji) WRITE(numout,*) 'Snow mass on ice :', sum_snow_nobio(ji) WRITE(numout,*) 'Melt water :', tot_melt(ji) WRITE(numout,*) 'evapwet : ', vevapwet(ji,:) WRITE(numout,*) 'transpir : ', transpir(ji,:) WRITE(numout,*) 'evapnu, evapsno : ', vevapnu(ji), vevapsno(ji) WRITE(numout,*) 'drainage : ', drainage(ji) STOP 'in hydrolc_waterbal' ENDIF ! ENDDO ! ! Transfer the total water amount at the end of the current timestep top the begining of the next one. ! tot_water_beg = tot_water_end ! END SUBROUTINE hydrolc_waterbal ! ! This routine computes the changes in soil moisture and interception storage for the ALMA outputs ! SUBROUTINE hydrolc_alma (kjpindex, index, first_call, qsintveg, snow, snow_nobio, soilwet) ! INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map LOGICAL, INTENT (in) :: first_call !! At which time is this routine called ? ! REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow water equivalent REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilwet !! Soil wetness REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Water balance on ice, lakes, .. [Kg/m^2] ! ! LOCAL ! INTEGER(i_std) :: ji REAL(r_std) :: watveg ! ! ! IF ( first_call ) THEN tot_watveg_beg(:) = zero tot_watsoil_beg(:) = zero snow_beg(:) = zero ! DO ji = 1, kjpindex watveg = SUM(qsintveg(ji,:)) tot_watveg_beg(ji) = watveg tot_watsoil_beg(ji) = mean_bqsb(ji) + mean_gqsb(ji) snow_beg(ji) = snow(ji)+ SUM(snow_nobio(ji,:)) ENDDO ! tot_watveg_end(:) = tot_watveg_beg(:) tot_watsoil_end(:) = tot_watsoil_beg(:) snow_end(:) = snow_beg(:) RETURN ENDIF ! ! Calculate the values for the end of the time step ! tot_watveg_end(:) = zero tot_watsoil_end(:) = zero snow_end(:) = zero delintercept(:) = zero delsoilmoist(:) = zero delswe(:) = zero ! DO ji = 1, kjpindex watveg = SUM(qsintveg(ji,:)) tot_watveg_end(ji) = watveg tot_watsoil_end(ji) = mean_bqsb(ji) + mean_gqsb(ji) snow_end(ji) = snow(ji)+ SUM(snow_nobio(ji,:)) ! delintercept(ji) = tot_watveg_end(ji) - tot_watveg_beg(ji) delsoilmoist(ji) = tot_watsoil_end(ji) - tot_watsoil_beg(ji) delswe(ji) = snow_end(ji) - snow_beg(ji) ! ! ENDDO ! ! ! Transfer the total water amount at the end of the current timestep top the begining of the next one. ! tot_watveg_beg = tot_watveg_end tot_watsoil_beg = tot_watsoil_end snow_beg(:) = snow_end(:) ! DO ji = 1,kjpindex soilwet(ji) = tot_watsoil_end(ji) / mx_eau_var(ji) ENDDO ! END SUBROUTINE hydrolc_alma !- SUBROUTINE hydrolc_hdiff(kjpindex, dtradia, veget, ruu_ch, gqsb, bqsb, dsg, dss, dsp) INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size REAL(r_std), INTENT (in) :: dtradia !! time step (s) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ruu_ch !! Quantite d'eau maximum REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: gqsb !! Hauteur d'eau dans le reservoir de surface REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: bqsb !! Hauteur d'eau dans le reservoir profond REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsg !! Hauteur du reservoir de surface REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dss !! Hauteur au dessus du reservoir de surface REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsp !! Hauteur au dessus du reservoir profond REAL(r_std), DIMENSION (kjpindex) :: bqsb_mean REAL(r_std), DIMENSION (kjpindex) :: gqsb_mean REAL(r_std), DIMENSION (kjpindex) :: dss_mean REAL(r_std), DIMENSION (kjpindex) :: vegtot REAL(r_std) :: x INTEGER(i_std) :: ji,jv REAL(r_std), SAVE :: tau_hdiff LOGICAL, SAVE :: firstcall=.TRUE. IF ( firstcall ) THEN !Config Key = HYDROL_TAU_HDIFF !Config Desc = time scale (s) for horizontal diffusion of water !Config Def = one_day !Config If = HYDROL_OK_HDIFF !Config Help = Defines how fast diffusion occurs horizontally between !Config the individual PFTs' water reservoirs. If infinite, no !Config diffusion. tau_hdiff = one_day CALL getin_p('HYDROL_TAU_HDIFF',tau_hdiff) WRITE (numout,*) 'Hydrol: Horizontal diffusion, tau (s)=',tau_hdiff firstcall = .FALSE. ENDIF ! Calculate mean values ! could be done with SUM instruction but this kills vectorization ! bqsb_mean(:) = zero gqsb_mean(:) = zero dss_mean(:) = zero vegtot(:) = zero ! DO jv = 1, nvm DO ji = 1, kjpindex bqsb_mean(ji) = bqsb_mean(ji) + veget(ji,jv)*bqsb(ji,jv) gqsb_mean(ji) = gqsb_mean(ji) + veget(ji,jv)*gqsb(ji,jv) dss_mean(ji) = dss_mean(ji) + veget(ji,jv)*dss(ji,jv) vegtot(ji) = vegtot(ji) + veget(ji,jv) ENDDO ENDDO DO ji = 1, kjpindex IF (vegtot(ji) .GT. zero) THEN bqsb_mean(ji) = bqsb_mean(ji)/vegtot(ji) gqsb_mean(ji) = gqsb_mean(ji)/vegtot(ji) dss_mean(ji) = dss_mean(ji)/vegtot(ji) ENDIF ENDDO ! relax values towards mean. ! x = MAX( zero, MIN( dtradia/tau_hdiff, un ) ) ! DO jv = 1, nvm DO ji = 1, kjpindex ! bqsb(ji,jv) = (un-x) * bqsb(ji,jv) + x * bqsb_mean(ji) gqsb(ji,jv) = (un-x) * gqsb(ji,jv) + x * gqsb_mean(ji) dss(ji,jv) = (un-x) * dss(ji,jv) + x * dss_mean(ji) ! IF (gqsb(ji,jv) .LT. min_sechiba) THEN dsg(ji,jv) = zero ELSE dsg(ji,jv) = (dss(ji,jv) * ruu_ch(ji) + gqsb(ji,jv)) / ruu_ch(ji) ENDIF dsp(ji,jv) = dpu_cste - bqsb(ji,jv) / ruu_ch(ji) ! ENDDO ENDDO END SUBROUTINE hydrolc_hdiff ! END MODULE hydrolc