!! !! This module computes hydrologic processes on continental points. !! !! @author Marie-Alice Foujols and Jan Polcher !! @Version : $Revision: 1.36 $, $Date: 2009/01/07 13:39:45 $ !! !! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/hydrol.f90,v 1.36 2009/01/07 13:39:45 ssipsl Exp $ !! IPSL (2006) !! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC !! MODULE hydrol ! ! ! routines called : restput, restget ! USE ioipsl ! ! modules used : USE constantes USE constantes_soil USE constantes_veg USE sechiba_io ! for debug : USE grid IMPLICIT NONE ! public routines : ! hydrol PRIVATE PUBLIC :: hydrol_main,hydrol_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=.TRUE. !! The check the water balance LOGICAL, SAVE :: check_cwrr=.TRUE. !! The check the water balance CHARACTER(LEN=80) , SAVE :: file_ext !! Extention for I/O filename CHARACTER(LEN=80) , SAVE :: var_name !! To store variables names for I/O REAL(r_std), PARAMETER :: drain_rest_cste = 15.0 !! time constant in days to return to free drainage after return flow REAL(r_std), PARAMETER :: allowed_err = 1.0E-8_r_std REAL(r_std), PARAMETER :: dmcs = 0.002 !! Allowed moisture above mcs (boundary conditions) REAL(r_std), PARAMETER :: dmcr = 0.002 !! Allowed moisture below mcr (boundary conditions) ! one dimension array allocated, computed, saved and got in hydrol module 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^Q ! array allocated, computed, saved and got in hydrol module INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mask_veget !! zero/one when veget fraction is zero/higher INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mask_soiltype !! zero/one where soil fraction is zero/higher INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mask_corr_veg_soil !! zero/one where veg frac on a soil type is zero/higher INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mask_return !! zero/one where there is no/is returnflow INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: index_nsat !! Indices of the non-saturated layers INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: index_sat !! Indices of the saturated layers INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: n_nsat !! Number of non-saturated layers INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: n_sat !! Number of saturated layers INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: nslme !! last efficient layer REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: humrelv !! humrel for each soil type REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: vegstressv !! vegstress for each soil type REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:,:):: us !! relative humidity REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: precisol !! Eau tombee sur le sol REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: precisol_ns !! Eau tombee sur le sol par type de sol REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ae_ns !! Evaporation du sol nu par type de sol REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: evap_bare_lim_ns !! limitation of bare soil evaporation on each soil column (used to deconvoluate vevapnu) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: free_drain_coef !! Coefficient for free drainage at bottom REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: rootsink !! stress racinaire par niveau et type de 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 (:) :: subsinksoil !! Excess of sublimation as a sink for the soil 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 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mx_eau_var ! arrays used by cwrr scheme REAL(r_std), SAVE, DIMENSION (nslm+1,nstm) :: zz !! REAL(r_std), SAVE, DIMENSION (nslm+1,nstm) :: dz !! REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: mc_lin !! REAL(r_std), SAVE, DIMENSION (nstm) :: v1r !! Residual soil water content of the first layer REAL(r_std), SAVE, DIMENSION (nstm) :: vBs !! Saturated soil water content of the bottom layer REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: k_lin !! REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: d_lin !! REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: a_lin !! REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: b_lin !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: humtot !! (:) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: flux !! (:) LOGICAL, ALLOCATABLE, SAVE, DIMENSION (:) :: resolv !! (:) !! linarization coefficients of hydraulic conductivity K REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: a !! (:,nslm) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: b !! !! linarization coefficients of hydraulic diffusivity D REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: d !! !! matrix coefficients REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: e !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: f !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: g1 !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ep !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: fp !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gp !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rhs !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: srhs !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gam !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc !! (:,nstm) Total moisture content (mm) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tmcs !! (nstm) Total moisture constent at saturation (mm) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter !! (:,nstm) Total moisture in the litter by soil type REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tmc_litt_mea !! Total moisture in the litter over the grid REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_wilt !! (:,nstm) Moisture of litter at wilt pt REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_field !! (:,nstm) Moisture of litter at field cap. REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_res !! (:,nstm) Moisture of litter at residual moisture. REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_sat !! (:,nstm) Moisture of litter at saturatiion REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_awet !! (:,nstm) Moisture of litter at mc_awet REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_adry !! (:,nstm) Moisture of litter at mc_dry REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tmc_litt_wet_mea !! Total moisture in the litter over the grid below which albedo is fixed REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tmc_litt_dry_mea !! Total moisture in the litter over the grid above which albedo is fixed REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: v1 !! (:) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vB !! (:) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: qflux00 !! flux at the top of the soil column !! par type de sol : REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ru_ns !! ruissellement REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dr_ns !! drainage REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tr_ns !! transpiration REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: corr_veg_soil !! (:,nvm,nstm) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: corr_veg_soil_max !! (:,nvm,nstm) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc !! (:,nslm,nstm) m³ x m³ REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soilmoist !! (:,nslm) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: soil_wet !! soil wetness REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soil_wet_litter !! soil wetness of the litter REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: qflux !! fluxes between the soil layers REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: tmat !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: stmat !! LOGICAL, SAVE :: interpol_diag=.FALSE. 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 hydrol_snow for snow process (including age of snow) !! - call hydrol_canop for canopy process !! - call hydrol_soil for bare soil process !! !! @call hydrol_snow !! @call hydrol_canop !! @call hydrol_soil !! SUBROUTINE hydrol_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, & & index, indexveg, indexsoil, indexlayer,& & temp_sol_new, runoff, 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, drysoil_frac, evapot, evapot_penm, evap_bare_lim, & & shumdiag, litterhumdiag, soilcap, soiltype, 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 for veg INTEGER(i_std),DIMENSION (kjpindex*nstm), INTENT (in):: indexsoil !! Indeces of the points on the 3D map for soil INTEGER(i_std),DIMENSION (kjpindex*nslm), INTENT (in):: indexlayer !! Indeces of the points on the 3D map for soil layers 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,nstm), INTENT (in) :: soiltype !! Map of soil types 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_penm !! Soil Potential Evaporation Correction ! modified fields REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: evap_bare_lim !! 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 ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: runoff !! 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), INTENT (out) :: drysoil_frac !! function of litter wetness 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) :: tot_melt !! Total melt REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintveg !! Water on vegetation due to interception ! ! local declaration ! INTEGER(i_std) :: jst, jsl REAL(r_std),DIMENSION (kjpindex) :: soilwet !! A temporary diagnostic of soil wetness REAL(r_std),DIMENSION (kjpindex) :: snowdepth !! Depth of snow layer ! ! do initialisation ! IF (l_first_hydrol) THEN IF (long_print) WRITE (numout,*) ' l_first_hydrol : call hydrol_init ' CALL hydrol_init (kjit, ldrestart_read, kjpindex, index, rest_id, veget, soiltype, humrel,& & vegstress, snow, snow_age, snow_nobio, snow_nobio_age, qsintveg) CALL hydrol_var_init (kjpindex, veget, soiltype, mx_eau_var, shumdiag, litterhumdiag, & & drysoil_frac, evap_bare_lim) ! ! If we check the water balance we first save the total amount of water ! IF (check_waterbal) THEN CALL hydrol_waterbal(kjpindex, index, .TRUE., dtradia, veget, & & totfrac_nobio, qsintveg, snow, snow_nobio,& & precip_rain, precip_snow, returnflow, irrigation, tot_melt, & & vevapwet, transpir, vevapnu, vevapsno, runoff,drainage) ENDIF ! IF (almaoutput) THEN CALL hydrol_alma(kjpindex, index, .TRUE., qsintveg, snow, snow_nobio, soilwet) ENDIF 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 ' DO jst=1,nstm ! var_name= "mc_1" ... "mc_3" WRITE (var_name,"('moistc_',i1)") jst CALL restput_p(rest_id, var_name, nbp_glo, nslm, 1, kjit, mc(:,:,jst), 'scatter', nbp_glo, index_g) END DO ! DO jst=1,nstm DO jsl=1,nslm ! var_name= "us_1_01" ... "us_3_11" WRITE (var_name,"('us_',i1,'_',i2.2)") jst,jsl CALL restput_p(rest_id, var_name, nbp_glo,nvm, 1,kjit,us(:,:,jst,jsl),'scatter',nbp_glo,index_g) END DO END DO ! var_name= 'free_drain_coef' CALL restput_p(rest_id, var_name, nbp_glo, nstm, 1, kjit, free_drain_coef, 'scatter', nbp_glo, index_g) ! var_name= 'ae_ns' CALL restput_p(rest_id, var_name, nbp_glo, nstm, 1, kjit, ae_ns, '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= '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 ! ! shared time step ! IF (long_print) WRITE (numout,*) 'hydrol pas de temps = ',dtradia ! ! computes snow ! CALL hydrol_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 hydrol_vegupd(kjpindex, veget, veget_max, soiltype, qsintveg,resdist) ! CALL hydrol_canop(kjpindex, precip_rain, vevapwet, veget, qsintmax, qsintveg,precisol,tot_melt) ! computes hydro_soil ! CALL hydrol_soil(kjpindex, dtradia, veget, veget_max, soiltype, transpir, vevapnu, evapot, & & evapot_penm, runoff, drainage, returnflow, irrigation, & & tot_melt,evap_bare_lim, shumdiag, litterhumdiag, humrel, vegstress, drysoil_frac) ! ! If we check the water balance we end with the comparison of total water change and fluxes ! IF (check_waterbal) THEN CALL hydrol_waterbal(kjpindex, index, .FALSE., dtradia, veget, totfrac_nobio, & & qsintveg, snow,snow_nobio, precip_rain, precip_snow, returnflow, & & irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno, runoff, drainage) ENDIF ! ! If we use the ALMA standards ! IF (almaoutput) THEN CALL hydrol_alma(kjpindex, index, .FALSE., qsintveg, snow, snow_nobio, soilwet) ENDIF ! IF ( .NOT. almaoutput ) THEN DO jst=1,nstm ! var_name= "mc_1" ... "mc_3" WRITE (var_name,"('moistc_',i1)") jst CALL histwrite(hist_id, trim(var_name), kjit,mc(:,:,jst), kjpindex*nslm, indexlayer) ! var_name= "vegetsoil_1" ... "vegetsoil_3" WRITE (var_name,"('vegetsoil_',i1)") jst CALL histwrite(hist_id, trim(var_name), kjit,corr_veg_soil(:,:,jst), kjpindex*nvm, indexveg) ENDDO CALL histwrite(hist_id, 'evapnu_soil', kjit, ae_ns, kjpindex*nstm, indexsoil) CALL histwrite(hist_id, 'drainage_soil', kjit, dr_ns, kjpindex*nstm, indexsoil) CALL histwrite(hist_id, 'transpir_soil', kjit, tr_ns, kjpindex*nstm, indexsoil) CALL histwrite(hist_id, 'runoff_soil', kjit, ru_ns, kjpindex*nstm, indexsoil) CALL histwrite(hist_id, 'humtot_soil', kjit, tmc, kjpindex*nstm, indexsoil) CALL histwrite(hist_id, 'humtot', kjit, humtot, kjpindex, index) CALL histwrite(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg) CALL histwrite(hist_id, 'drainage', kjit, drainage, kjpindex, index) CALL histwrite(hist_id, 'runoff', kjit, runoff, 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) IF ( hist2_id > 0 ) THEN DO jst=1,nstm ! var_name= "mc_1" ... "mc_3" WRITE (var_name,"('moistc_',i1)") jst CALL histwrite(hist2_id, trim(var_name), kjit,mc(:,:,jst), kjpindex*nslm, indexlayer) ! var_name= "vegetsoil_1" ... "vegetsoil_3" WRITE (var_name,"('vegetsoil_',i1)") jst CALL histwrite(hist2_id, trim(var_name), kjit,corr_veg_soil(:,:,jst), kjpindex*nvm, indexveg) ENDDO CALL histwrite(hist2_id, 'evapnu_soil', kjit, ae_ns, kjpindex*nstm, indexsoil) CALL histwrite(hist2_id, 'drainage_soil', kjit, dr_ns, kjpindex*nstm, indexsoil) CALL histwrite(hist2_id, 'transpir_soil', kjit, tr_ns, kjpindex*nstm, indexsoil) CALL histwrite(hist2_id, 'runoff_soil', kjit, ru_ns, kjpindex*nstm, indexsoil) CALL histwrite(hist2_id, 'humtot_soil', kjit, tmc, kjpindex*nstm, indexsoil) CALL histwrite(hist2_id, 'humtot', kjit, humtot, kjpindex, index) CALL histwrite(hist2_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg) CALL histwrite(hist2_id, 'drainage', kjit, drainage, kjpindex, index) CALL histwrite(hist2_id, 'runoff', kjit, runoff, 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) ENDIF 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, runoff, 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, soilmoist, kjpindex*nslm, indexlayer) CALL histwrite(hist_id, 'SoilWet', kjit, soilwet, kjpindex, index) ! 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) ! IF ( hist2_id > 0 ) THEN 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, runoff, 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, soilmoist, kjpindex*nslm, indexlayer) 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,*) ' hydrol_main Done ' END SUBROUTINE hydrol_main !! Algorithm: !! - dynamic allocation for local array !! - _restart_ file reading for HYDROLOGIC variables !! SUBROUTINE hydrol_init(kjit, ldrestart_read, kjpindex, index, rest_id, veget, soiltype, 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 REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types ! 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, ierror, ipdt INTEGER(i_std) :: ji, jv, jst, jsl, ik ! initialisation IF (l_first_hydrol) THEN l_first_hydrol=.FALSE. ELSE WRITE (numout,*) ' l_first_hydrol false . we stop ' STOP 'hydrol_init' ENDIF ! make dynamic allocation with good dimension ! one dimension array allocation with possible restart value ALLOCATE (mask_corr_veg_soil(kjpindex,nvm,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mask_corr_veg_soil allocation. We stop. We need kjpindex words = ',kjpindex*nvm*nstm STOP 'hydrol_init' END IF ALLOCATE (mask_veget(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mask_veget allocation. We stop. We need kjpindex words = ',kjpindex*nvm STOP 'hydrol_init' END IF ALLOCATE (mask_soiltype(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (mask_return(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (index_nsat(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (index_sat(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (n_nsat(nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',nstm STOP 'hydrol_init' END IF ALLOCATE (n_sat(nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',nstm STOP 'hydrol_init' END IF ALLOCATE (nslme(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (humrelv(kjpindex,nvm,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in humrelv allocation. We stop. We need kjpindex words = ',kjpindex*nvm*nstm STOP 'hydrol_init' END IF ALLOCATE (vegstressv(kjpindex,nvm,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in vegstressv allocation. We stop. We need kjpindex words = ',kjpindex*nvm*nstm STOP 'hydrol_init' END IF ALLOCATE (us(kjpindex,nvm,nstm,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in us allocation. We stop. We need kjpindex words = ',kjpindex*nvm*nstm*nslm STOP 'hydrol_init' END IF ALLOCATE (precisol(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindex*nvm STOP 'hydrol_init' END IF ALLOCATE (precisol_ns(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in precisol_ns allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (free_drain_coef(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in free_drain_coef allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (ae_ns(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in ae_ns allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (evap_bare_lim_ns(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in evap_bare_lim_ns allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (rootsink(kjpindex,nslm,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in rootsink allocation. We stop. We need kjpindex words = ',kjpindex*nslm*nstm STOP 'hydrol_init' END IF ALLOCATE (subsnowveg(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in subsnowveg allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (subsnownobio(kjpindex,nnobio),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in subsnownobio allocation. We stop. We need kjpindex words = ',kjpindex*nnobio STOP 'hydrol_init' END IF ALLOCATE (snowmelt(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in snowmelt allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (icemelt(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in icemelt allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (subsinksoil(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in subsinksoil allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF 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 'hydrol_init' END IF ALLOCATE (vegtot(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in vegtot allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF 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 'hydrol_init' END IF ALLOCATE (humtot(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in humtot allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (flux(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in flux allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (resolv(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in resolv allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (a(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in a allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (b(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in b allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (d(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in d allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (e(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in e allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (f(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in f allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (g1(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in g1 allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (ep(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in ep allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (fp(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in fp allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (gp(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in gp allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (rhs(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in rhs allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (srhs(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in srhs allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (gam(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in gam allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (tmc(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (tmcs(nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmcs allocation. We stop. We need kjpindex words = ',nstm STOP 'hydrol_init' END IF ALLOCATE (tmc_litter(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litter allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (tmc_litt_mea(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litt_mea allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (tmc_litter_res(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litter_res allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (tmc_litter_wilt(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litter_wilt allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (tmc_litter_field(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litter_field allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (tmc_litter_sat(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litter_sat allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (tmc_litter_awet(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litter_awet allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (tmc_litter_adry(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litter_adry allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (tmc_litt_wet_mea(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litt_wet_mea allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (tmc_litt_dry_mea(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmc_litt_dry_mea allocation. We stop. We need kjpindex words = ',kjpindex STOP 'hydrol_init' END IF ALLOCATE (v1(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in v1 allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (vB(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in vB allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (qflux00(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in qflux00 allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (ru_ns(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in ru_ns allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ru_ns(:,:) = zero ALLOCATE (dr_ns(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in dr_ns allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF dr_ns(:,:) = zero ALLOCATE (tr_ns(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tr_ns allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (corr_veg_soil(kjpindex,nvm,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in corr_veg_soil allocation. We stop. We need kjpindex words = ',kjpindex*nvm*nstm STOP 'hydrol_init' END IF ALLOCATE (corr_veg_soil_max(kjpindex,nvm,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in corr_veg_soil_max allocation. We stop. We need kjpindex words = ',kjpindex*nvm*nstm STOP 'hydrol_init' END IF ALLOCATE (mc(kjpindex,nslm,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in mc allocation. We stop. We need kjpindex words = ',kjpindex*nslm*nstm STOP 'hydrol_init' END IF ALLOCATE (soilmoist(kjpindex,nslm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in soilmoist allocation. We stop. We need kjpindex words = ',kjpindex*nslm STOP 'hydrol_init' END IF ALLOCATE (soil_wet(kjpindex,nslm,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in soil_wet allocation. We stop. We need kjpindex words = ',kjpindex*nslm*nstm STOP 'hydrol_init' END IF ALLOCATE (soil_wet_litter(kjpindex,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in soil_wet allocation. We stop. We need kjpindex words = ',kjpindex*nstm STOP 'hydrol_init' END IF ALLOCATE (qflux(kjpindex,nslm,nstm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in qflux allocation. We stop. We need kjpindex words = ',kjpindex*nslm*nstm STOP 'hydrol_init' END IF ALLOCATE (tmat(kjpindex,nslm,3),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in tmat allocation. We stop. We need kjpindex words = ',kjpindex*nslm*trois STOP 'hydrol_init' END IF ALLOCATE (stmat(kjpindex,nslm,3),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in stmat allocation. We stop. We need kjpindex words = ',kjpindex*nslm*trois STOP 'hydrol_init' END IF ! ! 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 'hydrol_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 'hydrol_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 'hydrol_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 'hydrol_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 'hydrol_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 'hydrol_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 'hydrol_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 'hydrol_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 'hydrol_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 'hydrol_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 'hydrol_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' IF (is_root_prc) CALL ioconf_setatt('UNITS', '-') ! DO jst=1,nstm ! var_name= "mc_1" ... "mc_3" WRITE (var_name,"('moistc_',I1)") jst CALL ioconf_setatt('LONG_NAME',var_name) CALL restget_p (rest_id, var_name, nbp_glo, nslm , 1, kjit, .TRUE., mc(:,:,jst), "gather", nbp_glo, index_g) END DO ! CALL ioconf_setatt('UNITS', '-') DO jst=1,nstm DO jsl=1,nslm ! var_name= "us_1_01" ... "us_3_11" WRITE (var_name,"('us_',i1,'_',i2.2)") jst,jsl CALL ioconf_setatt('LONG_NAME',var_name) CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., us(:,:,jst,jsl), "gather", nbp_glo, index_g) END DO END DO ! var_name= 'free_drain_coef' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Coefficient for free drainage at bottom of soil') CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., free_drain_coef, "gather", nbp_glo, index_g) ! var_name= 'ae_ns' CALL ioconf_setatt('UNITS', 'kg/m^2') CALL ioconf_setatt('LONG_NAME','Bare soil evap on each soil type') CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., ae_ns, "gather", nbp_glo, index_g) ! 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) ! 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) ! 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= '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= '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_MOISTURE_CONTENT !Config Desc = Soil moisture on each soil tile and levels !Config Def = 0.3 !Config Help = The initial value of mc 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 (mc, val_exp, 'HYDROL_MOISTURE_CONTENT', 0.3_r_std) ! !Config Key = US_INIT !Config Desc = US_NVM_NSTM_NSLM !Config Def = 0.0 !Config Help = The initial value of us 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. ! DO jsl=1,nslm CALL setvar_p (us(:,:,:,jsl), val_exp, 'US_INIT', 0.0_r_std) ENDDO ! !Config Key = FREE_DRAIN_COEF !Config Desc = Coefficient for free drainage at bottom !Config Def = 1.0, 1.0, 1.0 !Config Help = The initial value of free drainage 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 (free_drain_coef, val_exp, 'FREE_DRAIN_COEF', free_drain_max) ! !Config Key = EVAPNU_SOIL !Config Desc = Bare soil evap on each soil if not found in restart !Config Def = 0.0 !Config Help = The initial value of bare soils evap 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 (ae_ns, val_exp, 'EVAPNU_SOIL', 0.0_r_std) ! !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', 0.0_r_std) ! !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', 0.0_r_std) ! !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', 0.0_r_std) ! !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', 0.0_r_std) ! ! ! !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', 0.0_r_std) ! ! 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 ! ! ! compute the masks for veget ! mask_veget(:,:) = MIN( un, MAX(zero,veget(:,:))) ! mask_soiltype(:,:) = MIN( un, MAX(zero,soiltype(:,:))) ! mask_corr_veg_soil(:,:,:) = MIN( un, MAX(zero,corr_veg_soil(:,:,:))) mask_veget(:,:) = 0 mask_soiltype(:,:) = 0 mask_corr_veg_soil(:,:,:) = 0 mask_return(:) = 0 index_nsat(:,:) = 0 index_sat(:,:) = 0 n_nsat(:) = 1 n_sat(:) = 0 nslme(:,:) = nslm DO ji = 1, kjpindex DO jst = 1, nstm IF(soiltype(ji,jst) .GT. min_sechiba) THEN mask_soiltype(ji,jst) = 1 ENDIF END DO DO jv = 1, nvm IF(veget(ji,jv) .GT. min_sechiba) THEN mask_veget(ji,jv) = 1 ENDIF DO jst = 1, nstm IF(corr_veg_soil(ji,jv,jst) .GT. min_sechiba) THEN mask_corr_veg_soil(ji,jv,jst) = 1 ENDIF END DO END DO ! WRITE(numout,*) 'mask: soiltype,mask_soiltype',soiltype(ji,:),mask_soiltype(ji,:) END DO ! set humrelv from us humrelv(:,:,:) = SUM(us,dim=4) vegstressv(:,:,:) = SUM(us,dim=4) ! set humrel from humrelv humrel(:,:) = zero DO jst=1,nstm DO jv=1,nvm DO ji=1,kjpindex vegstress(ji,jv)=vegstress(ji,jv) + vegstressv(ji,jv,jst) * soiltype(ji,jst) ! IF(veget(ji,jv).NE.zero) THEN humrel(ji,jv)=humrel(ji,jv) + humrelv(ji,jv,jst) * & & soiltype(ji,jst) humrel(ji,jv)=MAX(humrel(ji,jv), zero)* mask_veget(ji,jv) ! ELSE ! humrel(ji,jv)= zero ! ENDIF END DO END DO END DO ! vegstress(:,:)=humrel(:,:) ENDIF ! ! IF (long_print) WRITE (numout,*) ' hydrol_init done ' ! END SUBROUTINE hydrol_init ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE hydrol_clear() l_first_hydrol=.TRUE. IF (ALLOCATED (mask_veget)) DEALLOCATE (mask_veget) IF (ALLOCATED (mask_soiltype)) DEALLOCATE (mask_soiltype) IF (ALLOCATED (mask_corr_veg_soil)) DEALLOCATE (mask_corr_veg_soil) IF (ALLOCATED (mask_return)) DEALLOCATE (mask_return) IF (ALLOCATED (index_nsat)) DEALLOCATE (index_nsat) IF (ALLOCATED (index_sat)) DEALLOCATE (index_sat) IF (ALLOCATED (n_nsat)) DEALLOCATE (n_nsat) IF (ALLOCATED (n_sat)) DEALLOCATE (n_sat) IF (ALLOCATED (nslme)) DEALLOCATE (nslme) IF (ALLOCATED (humrelv)) DEALLOCATE (humrelv) IF (ALLOCATED (vegstressv)) DEALLOCATE (vegstressv) IF (ALLOCATED (us)) DEALLOCATE (us) IF (ALLOCATED (precisol)) DEALLOCATE (precisol) IF (ALLOCATED (precisol_ns)) DEALLOCATE (precisol_ns) IF (ALLOCATED (free_drain_coef)) DEALLOCATE (free_drain_coef) IF (ALLOCATED (ae_ns)) DEALLOCATE (ae_ns) IF (ALLOCATED (evap_bare_lim_ns)) DEALLOCATE (evap_bare_lim_ns) IF (ALLOCATED (rootsink)) DEALLOCATE (rootsink) IF (ALLOCATED (subsnowveg)) DEALLOCATE (subsnowveg) IF (ALLOCATED (subsnownobio)) DEALLOCATE (subsnownobio) IF (ALLOCATED (snowmelt)) DEALLOCATE (snowmelt) IF (ALLOCATED (icemelt)) DEALLOCATE (icemelt) IF (ALLOCATED (subsinksoil)) DEALLOCATE (subsinksoil) IF (ALLOCATED (mx_eau_var)) DEALLOCATE (mx_eau_var) IF (ALLOCATED (vegtot)) DEALLOCATE (vegtot) IF (ALLOCATED (resdist)) DEALLOCATE (resdist) 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) ! more allocation for cwrr scheme IF (ALLOCATED (v1)) DEALLOCATE (v1) IF (ALLOCATED (vB)) DEALLOCATE (vB) IF (ALLOCATED (humtot)) DEALLOCATE (humtot) IF (ALLOCATED (flux)) DEALLOCATE (flux) IF (ALLOCATED (resolv)) DEALLOCATE (resolv) IF (ALLOCATED (a)) DEALLOCATE (a) IF (ALLOCATED (b)) DEALLOCATE (b) IF (ALLOCATED (d)) DEALLOCATE (d) IF (ALLOCATED (e)) DEALLOCATE (e) IF (ALLOCATED (f)) DEALLOCATE (f) IF (ALLOCATED (g1)) DEALLOCATE (g1) IF (ALLOCATED (ep)) DEALLOCATE (ep) IF (ALLOCATED (fp)) DEALLOCATE (fp) IF (ALLOCATED (gp)) DEALLOCATE (gp) IF (ALLOCATED (rhs)) DEALLOCATE (rhs) IF (ALLOCATED (srhs)) DEALLOCATE (srhs) IF (ALLOCATED (gam)) DEALLOCATE (gam) IF (ALLOCATED (tmc)) DEALLOCATE (tmc) IF (ALLOCATED (tmcs)) DEALLOCATE (tmcs) IF (ALLOCATED (tmc_litter)) DEALLOCATE (tmc_litter) IF (ALLOCATED (tmc_litt_mea)) DEALLOCATE (tmc_litt_mea) IF (ALLOCATED (tmc_litter_res)) DEALLOCATE (tmc_litter_res) IF (ALLOCATED (tmc_litter_wilt)) DEALLOCATE (tmc_litter_wilt) IF (ALLOCATED (tmc_litter_field)) DEALLOCATE (tmc_litter_field) IF (ALLOCATED (tmc_litter_sat)) DEALLOCATE (tmc_litter_sat) IF (ALLOCATED (tmc_litter_awet)) DEALLOCATE (tmc_litter_awet) IF (ALLOCATED (tmc_litter_adry)) DEALLOCATE (tmc_litter_adry) IF (ALLOCATED (tmc_litt_wet_mea)) DEALLOCATE (tmc_litt_wet_mea) IF (ALLOCATED (tmc_litt_dry_mea)) DEALLOCATE (tmc_litt_dry_mea) IF (ALLOCATED (qflux00)) DEALLOCATE (qflux00) IF (ALLOCATED (ru_ns)) DEALLOCATE (ru_ns) IF (ALLOCATED (dr_ns)) DEALLOCATE (dr_ns) IF (ALLOCATED (tr_ns)) DEALLOCATE (tr_ns) IF (ALLOCATED (corr_veg_soil)) DEALLOCATE (corr_veg_soil) IF (ALLOCATED (corr_veg_soil_max)) DEALLOCATE (corr_veg_soil_max) IF (ALLOCATED (mc)) DEALLOCATE (mc) IF (ALLOCATED (soilmoist)) DEALLOCATE (soilmoist) IF (ALLOCATED (soil_wet)) DEALLOCATE (soil_wet) IF (ALLOCATED (soil_wet_litter)) DEALLOCATE (soil_wet_litter) IF (ALLOCATED (qflux)) DEALLOCATE (qflux) IF (ALLOCATED (tmat)) DEALLOCATE (tmat) IF (ALLOCATED (stmat)) DEALLOCATE (stmat) RETURN END SUBROUTINE hydrol_clear !! This routine initializes HYDROLOGIC variables !! - mx_eau_var SUBROUTINE hydrol_var_init (kjpindex, veget, soiltype, mx_eau_var, shumdiag, litterhumdiag, drysoil_frac, evap_bare_lim) ! 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,nstm), INTENT (in) :: soiltype !! Map of soil types ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: mx_eau_var !! 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) :: drysoil_frac !! function of litter humidity REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: evap_bare_lim !! ! local declaration REAL(r_std), DIMENSION (kjpindex) :: dpu_mean !! mean soil depth INTEGER(i_std) :: ji, jv, jd, jst, jsl REAL(r_std) :: m, frac ! ! initialisation mx_eau_var(:) = zero dpu_mean(:)= zero ! DO ji = 1,kjpindex DO jst = 1,nstm dpu_mean(ji)=dpu_mean(ji)+dpu(jst)*soiltype(ji,jst) END DO END DO ! DO ji = 1,kjpindex DO jst = 1,nstm mx_eau_var(ji) = mx_eau_var(ji) + soiltype(ji,jst)*& & dpu(jst)*mille*mcs(jst) END DO END DO DO ji = 1,kjpindex IF (vegtot(ji) .LE. zero) THEN mx_eau_var(ji) = mx_eau_eau*deux ENDIF END DO ! ! Calcul the matrix coef for dublin model: ! pice-wise linearised hydraulic conductivity k_lin=alin * mc_lin + b_lin ! and diffusivity d_lin in each interval of mc, called mc_lin, ! between imin, for residual mcr, ! and imax for saturation mcs. ! DO jst=1,nstm m = un - un / nvan(jst) ! WRITE(numout,*) 'jst',jst,imin,imax mc_lin(imin,jst)=mcr(jst) mc_lin(imax,jst)=mcs(jst) tmcs(jst)=dpu(jst)* mille*mcs(jst) zz(1,jst) = zero dz(1,jst) = zero DO jsl=2,nslm zz(jsl,jst) = dpu(jst)* mille*((2**(jsl-1))-1)/ ((2**(nslm-1))-1) dz(jsl,jst) = zz(jsl,jst)-zz(jsl-1,jst) ! WRITE(numout,*) 'jsl, zz,dz',jsl, dpu(jst),zz(jsl,jst),dz(jsl,jst) ENDDO zz(nslm+1,jst) = zz(nslm,jst) dz(nslm+1,jst) = zero DO ji= imin+1, imax-1 mc_lin(ji,jst) = mcr(jst) + (ji-imin)*(mcs(jst)-mcr(jst))/(imax-imin) ENDDO DO ji = imin,imax-1 frac=MIN(un,(mc_lin(ji,jst)-mcr(jst))/(mcs(jst)-mcr(jst))) k_lin(ji,jst) = ks(jst) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2 frac=MIN(un,(mc_lin(ji+1,jst)-mcr(jst))/(mcs(jst)-mcr(jst))) k_lin(ji+1,jst) = ks(jst) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2 a_lin(ji,jst) = (k_lin(ji+1,jst)-k_lin(ji,jst)) / (mc_lin(ji+1,jst)-mc_lin(ji,jst)) b_lin(ji,jst) = k_lin(ji,jst) - a_lin(ji,jst)*mc_lin(ji,jst) !- Il faudrait ici definir a et b pour mc > mcs, et mc < mcr car c'est un cas auquel on peut etre confronte... a reflechir IF(ji.NE.imin.AND.ji.NE.imax-1) THEN frac=MIN(un,(mc_lin(ji,jst)-mcr(jst))/(mcs(jst)-mcr(jst))) d_lin(ji,jst) =(k_lin(ji,jst) / (avan(jst)*m*nvan(jst))) * & & ( (frac**(-un/m))/(mc_lin(ji,jst)-mcr(jst)) ) * & & ( frac**(-un/m) -un ) ** (-m) frac=MIN(un,(mc_lin(ji+1,jst)-mcr(jst))/(mcs(jst)-mcr(jst))) d_lin(ji+1,jst) =(k_lin(ji+1,jst) / (avan(jst)*m*nvan(jst)))*& & ( (frac**(-un/m))/(mc_lin(ji+1,jst)-mcr(jst)) ) * & & ( frac**(-un/m) -un ) ** (-m) d_lin(ji,jst) = undemi * (d_lin(ji,jst)+d_lin(ji+1,jst)) ELSEIF(ji.EQ.imin) THEN d_lin(ji,jst) = zero ELSEIF(ji.EQ.imax-1) THEN frac=MIN(un,(mc_lin(ji,jst)-mcr(jst))/(mcs(jst)-mcr(jst))) d_lin(ji,jst) =(k_lin(ji,jst) / (avan(jst)*m*nvan(jst))) * & & ( (frac**(-un/m))/(mc_lin(ji,jst)-mcr(jst)) ) * & & ( frac**(-un/m) -un ) ** (-m) ENDIF ENDDO ENDDO ! Compute the litter humidity, shumdiag and fry litterhumdiag(:) = zero tmc_litt_mea(:) = zero tmc_litt_wet_mea(:) = zero tmc_litt_dry_mea(:) = zero shumdiag(:,:) = zero soilmoist(:,:) = zero humtot(:) = zero tmc(:,:) = zero ! Loop on soil types to compute the variables (ji,jst) DO jst=1,nstm ! the residual 1st layer soil moisture: v1r(jst) = dz(2,jst)*mcr(jst)/deux ! the saturated Bottom layer soil moisture: ! v A CALCULER SUR TOUT LE PROFIL (vs et vr egalement) vBs(jst) = dz(nslm,jst)*mcs(jst)/deux ! The total soil moisture for each soil type: DO ji=1,kjpindex tmc(ji,jst)= dz(2,jst) * ( trois*mc(ji,1,jst)+ mc(ji,2,jst))/huit END DO DO jsl=2,nslm-1 DO ji=1,kjpindex tmc(ji,jst) = tmc(ji,jst) + dz(jsl,jst) * ( trois*mc(ji,jsl,jst) + mc(ji,jsl-1,jst))/huit & & + dz(jsl+1,jst)*(trois*mc(ji,jsl,jst) + mc(ji,jsl+1,jst))/huit END DO END DO DO ji=1,kjpindex tmc(ji,jst) = tmc(ji,jst) + dz(nslm,jst) * (trois * mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit END DO ! The litter variables: DO ji=1,kjpindex tmc_litter(ji,jst) = dz(2,jst) * (trois*mc(ji,1,jst)+mc(ji,2,jst))/huit tmc_litter_wilt(ji,jst) = dz(2,jst) * mcw(jst) / deux tmc_litter_res(ji,jst) = dz(2,jst) * mcr(jst) / deux tmc_litter_field(ji,jst) = dz(2,jst) * mcf(jst) / deux tmc_litter_sat(ji,jst) = dz(2,jst) * mcs(jst) / deux tmc_litter_awet(ji,jst) = dz(2,jst) * mc_awet(jst) / deux tmc_litter_adry(ji,jst) = dz(2,jst) * mc_adry(jst) / deux END DO ! sum from level 1 to 4 DO jsl=2,4 DO ji=1,kjpindex tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl,jst) * & & ( trois*mc(ji,jsl,jst) + mc(ji,jsl-1,jst))/huit & & + dz(jsl+1,jst)*(trois*mc(ji,jsl,jst) + mc(ji,jsl+1,jst))/huit tmc_litter_wilt(ji,jst) = tmc_litter_wilt(ji,jst) + & &(dz(jsl,jst)+ dz(jsl+1,jst))*& & mcw(jst)/deux tmc_litter_res(ji,jst) = tmc_litter_res(ji,jst) + & &(dz(jsl,jst)+ dz(jsl+1,jst))*& & mcr(jst)/deux tmc_litter_sat(ji,jst) = tmc_litter_sat(ji,jst) + & &(dz(jsl,jst)+ dz(jsl+1,jst))* & & mcs(jst)/deux tmc_litter_field(ji,jst) = tmc_litter_field(ji,jst) + & & (dz(jsl,jst)+ dz(jsl+1,jst))* & & mcf(jst)/deux tmc_litter_awet(ji,jst) = tmc_litter_awet(ji,jst) + & &(dz(jsl,jst)+ dz(jsl+1,jst))* & & mc_awet(jst)/deux tmc_litter_adry(ji,jst) = tmc_litter_adry(ji,jst) + & & (dz(jsl,jst)+ dz(jsl+1,jst))* & & mc_adry(jst)/deux END DO END DO ! subsequent calcul of soil_wet_litter (tmc-tmcw)/(tmcf-tmcw) DO ji=1,kjpindex soil_wet_litter(ji,jst)=MIN(un, MAX(zero,& &(tmc_litter(ji,jst)-tmc_litter_wilt(ji,jst))/& & (tmc_litter_field(ji,jst)-tmc_litter_wilt(ji,jst)) )) END DO ! Soil wetness profiles (mc-mcw)/(mcs-mcw) DO ji=1,kjpindex soil_wet(ji,1,jst) = MIN(un, MAX(zero,& &(trois*mc(ji,1,jst) + mc(ji,2,jst) - quatre*mcw(jst))& & /(quatre*(mcs(jst)-mcw(jst))) )) humrelv(ji,1,jst) = zero END DO DO jsl=2,nslm-1 DO ji=1,kjpindex soil_wet(ji,jsl,jst) = MIN(un, MAX(zero,& & (trois*mc(ji,jsl,jst) + & & mc(ji,jsl-1,jst) *(dz(jsl,jst)/(dz(jsl,jst)+dz(jsl+1,jst))) & & + mc(ji,jsl+1,jst)*(dz(jsl+1,jst)/(dz(jsl,jst)+dz(jsl+1,jst))) & & - quatre*mcw(jst)) / (quatre*(mcs(jst)-mcw(jst))) )) END DO END DO DO ji=1,kjpindex soil_wet(ji,nslm,jst) = MIN(un, MAX(zero,& & (trois*mc(ji,nslm,jst) & & + mc(ji,nslm-1,jst)-quatre*mcw(jst))/(quatre*(mcs(jst)-mcw(jst))) )) END DO END DO ! loop on soil type !Now we compute the grid averaged values: DO jst=1,nstm DO ji=1,kjpindex humtot(ji) = humtot(ji) + soiltype(ji,jst) * tmc(ji,jst) litterhumdiag(ji) = litterhumdiag(ji) + & & soil_wet_litter(ji,jst) * soiltype(ji,jst) tmc_litt_wet_mea(ji) = tmc_litt_wet_mea(ji) + & & tmc_litter_awet(ji,jst)* soiltype(ji,jst) tmc_litt_dry_mea(ji) = tmc_litt_dry_mea(ji) + & & tmc_litter_adry(ji,jst) * soiltype(ji,jst) tmc_litt_mea(ji) = tmc_litt_mea(ji) + & & tmc_litter(ji,jst) * soiltype(ji,jst) END DO DO jsl=1,nbdl DO ji=1,kjpindex shumdiag(ji,jsl)= shumdiag(ji,jsl) + soil_wet(ji,jsl,jst) * & & ((mcs(jst)-mcw(jst))/(mcf(jst)-mcw(jst))) * & & soiltype(ji,jst) soilmoist(ji,jsl) = soilmoist(ji,jsl) + mc(ji,jsl,jst)*soiltype(ji,jst) shumdiag(ji,jsl) = MAX(MIN(shumdiag(ji,jsl), un), zero) END DO END DO END DO ! loop on soiltype ! ! ! DO ji=1,kjpindex drysoil_frac(ji) = un + MAX( MIN( (tmc_litt_dry_mea(ji) - tmc_litt_mea(ji)) / & & (tmc_litt_wet_mea(ji) - tmc_litt_dry_mea(ji)), zero), - un) END DO evap_bare_lim = zero !!$ IF ( COUNT(diaglev .EQ. undef_sechiba) > 0 ) THEN !!$ !!$ DO jsl=1,nbdl-1 !!$ diaglev(jsl) = zz(jsl,1) + dz(jsl+1,1)/deux !!$ END DO !!$ diaglev(nbdl) = zz(nbdl,1) !!$ interpol_diag = .FALSE. !!$ !!$ ENDIF IF (long_print) WRITE (numout,*) ' hydrol_var_init done ' END SUBROUTINE hydrol_var_init !! This routine computes snow processes !! SUBROUTINE hydrol_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 ! ! ! 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 subsinksoil(ji) = zero tot_melt(ji) = zero ENDDO ! ! 1. On vegetation ! DO ji=1,kjpindex ! ! 1.1. It is snowing ! snow(ji) = snow(ji) + (un - totfrac_nobio(ji))*precip_snow(ji) ! ! ! 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 ! ! ! 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 soil evaporation ! vevapnu(ji) = vevapnu(ji) + (subsnowveg(ji) - snow(ji)) IF( (un - totfrac_nobio(ji)).GT.min_sechiba) THEN subsinksoil (ji) = (subsnowveg(ji) - snow(ji))/ (un - totfrac_nobio(ji)) END IF ! 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 ! ! 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) = (1. - 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 WRITE(numout,*) 'hydrol_snow: WARNING! snow was negative and was reset to zero. ' ! END IF ENDIF ! ! 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' WRITE(numout,*) 'CANNOT TREAT SNOW ON THESE SURFACE TYPES' STOP 'in hydrol_snow' 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. min_sechiba ) 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,*) ' hydrol_snow done ' END SUBROUTINE hydrol_snow !! This routine computes canopy processes !! SUBROUTINE hydrol_canop (kjpindex, precip_rain, vevapwet, veget, qsintmax, & & qsintveg,precisol,tot_melt) ! ! 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 REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: tot_melt !! Total melt ! modified fields REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: qsintveg !! Water on vegetation due to interception REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: precisol !! Eau tombee sur le sol ! output fields ! ! local declaration ! INTEGER(i_std) :: ji, jv REAL(r_std), DIMENSION (kjpindex,nvm) :: zqsintvegnew LOGICAL, SAVE :: firstcall=.TRUE. REAL(r_std), SAVE, DIMENSION(nvm) :: throughfall_by_pft IF ( firstcall ) THEN !Config Key = PERCENT_THROUGHFALL_PFT !Config Desc = Percent by PFT of precip that is not intercepted by the canopy !Config Def = 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. !Config Help = During one rainfall event, PERCENT_THROUGHFALL_PFT% of the incident rainfall !Config will get directly to the ground without being intercepted, for each PFT. throughfall_by_pft = (/ 30., 30., 30., 30., 30., 30., 30., 30., 30., 30., 30., 30., 30. /) CALL getin_p('PERCENT_THROUGHFALL_PFT',throughfall_by_pft) throughfall_by_pft = throughfall_by_pft / 100. firstcall=.FALSE. ENDIF ! 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 END DO ! DO jv=1,nvm DO ji = 1, kjpindex IF (vegtot(ji).GT.min_sechiba) THEN precisol(ji,jv) = precisol(ji,jv)+tot_melt(ji)*veget(ji,jv)/vegtot(ji) ENDIF ENDDO END DO ! ! ! 1.4 swap qsintveg to the new value ! DO jv=1,nvm qsintveg(:,jv) = zqsintvegnew (:,jv) END DO IF (long_print) WRITE (numout,*) ' hydrol_canop done ' END SUBROUTINE hydrol_canop !! !! !! SUBROUTINE hydrol_vegupd(kjpindex, veget, veget_max, soiltype,qsintveg,resdist) ! ! The vegetation cover has changed and we need to adapt the reservoir distribution ! and the distribution of plants on different soil types. ! 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,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types : proportion of each soil type ! modified fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintveg !! Water on vegetation REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(inout) :: resdist !! Old vegetation map ! ! local declaration ! INTEGER(i_std) :: ji,jv,jst,jst_pref REAL(r_std), DIMENSION (kjpindex,nstm) :: soil_exist,soil_exist_max REAL(r_std), DIMENSION (kjpindex,nvm) :: veget_exist,veget_exist_max REAL(r_std), DIMENSION (kjpindex,nvm) :: qsintveg2 !! Water on vegetation due to interception over old veget REAL(r_std), DIMENSION (kjpindex,nvm) :: vmr !! variation of veget REAL(r_std), DIMENSION (kjpindex,nvm) :: qsdq REAL(r_std), DIMENSION(kjpindex) :: vegchtot,vtr, qstr, fra REAL(r_std), PARAMETER :: EPS1 = EPSILON(un) ! DO jv = 1, nvm DO ji = 1, kjpindex !mask ! vmr(ji,jv) = MAX ( EPSILON(un), MIN ( veget(ji,jv)-resdist(ji,jv) , MAX( EPSILON(un), veget(ji,jv)-resdist(ji,jv)) ) ) ! vmr(ji,jv) = MAX ( EPSILON(un), MAX( EPSILON(un), veget(ji,jv)-resdist(ji,jv)) ) ! IF(ABS(veget(ji,jv)-resdist(ji,jv)).gt.epsilon(un)) then ! WRITE(numout,*) '-----------------------------------------------' ! WRITE(numout,*) 'vmr,epsilon(un),veget,resdist',vmr(ji,jv),epsilon(un) ! WRITE(numout,*),veget(ji,jv),resdist(ji,jv) ! WRITE(numout,*) 'ABS(veget -resdist',ABS(veget(ji,jv)-resdist(ji,jv)) ! endif 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. min_sechiba) 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. min_sechiba ) THEN qsdq(ji,jv) = ABS(vmr(ji,jv)) * qsintveg2(ji,jv) ENDIF ENDDO ENDDO ! ! calculate water mass that we have to redistribute ! qstr(:) = zero vtr(:) = zero ! ! DO jv = 1, nvm DO ji = 1, kjpindex IF ( ( vegchtot(ji) .GT. min_sechiba ) .AND. ( vmr(ji,jv) .LT. -min_sechiba ) ) THEN 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. min_sechiba .AND. ABS(vtr(ji)) .GT. EPSILON(un)) THEN fra(ji) = vmr(ji,jv) / vtr(ji) IF ( vmr(ji,jv) .GT. min_sechiba) THEN qsintveg(ji,jv) = qsintveg(ji,jv) + fra(ji)* qstr(ji) ELSE qsintveg(ji,jv) = qsintveg(ji,jv) - qsdq(ji,jv) ENDIF ENDIF ENDDO ENDDO !MM underflow : DO jv = 1, nvm DO ji = 1, kjpindex IF ( ABS(qsintveg(ji,jv)) > 0. .AND. ABS(qsintveg(ji,jv)) < EPS1 ) THEN qsintveg(ji,jv) = EPS1 ENDIF ENDDO ENDDO ! Now that the work is done resdist needs an update ! DO jv = 1, nvm DO ji = 1, kjpindex resdist(ji,jv) = veget(ji,jv) ENDDO ENDDO ! Distribution of the vegetation depending on the soil type ! DO jst = 1, nstm ! ! DO ji = 1, kjpindex ! ! ! soil_exist(ji,jst)=zero ! IF (soiltype(ji,jst) .NE. zero) THEN ! soil_exist(ji,jst)=un ! soil_exist_max(ji,jst)=un ! ENDIF ! ! ENDDO ! ! ENDDO soil_exist(:,:) = mask_soiltype(:,:) soil_exist_max(:,:) = mask_soiltype(:,:) veget_exist(:,:) = zero veget_exist_max(:,:) = zero DO jv = 1, nvm DO ji = 1, kjpindex IF(vegtot(ji).GT.min_sechiba) THEN veget_exist(ji,jv)= veget(ji,jv)/vegtot(ji) veget_exist_max(ji,jv)= veget_max(ji,jv)/vegtot(ji) ENDIF ENDDO ENDDO ! Compute corr_veg_soil corr_veg_soil(:,:,:) = zero corr_veg_soil_max(:,:,:) = zero IF ( COUNT(pref_soil_veg .EQ. 0) > 0 ) THEN DO jst = 1, nstm DO jv = nvm, 1, -1 DO ji=1,kjpindex IF(vegtot(ji).GT.min_sechiba.AND.soiltype(ji,jst).GT.min_sechiba) THEN corr_veg_soil(ji,jv,jst) = veget(ji,jv)/vegtot(ji) corr_veg_soil_max(ji,jv,jst) = veget_max(ji,jv)/vegtot(ji) END IF END DO END DO END DO ELSE DO jst = 1, nstm DO jv = nvm, 1, -1 jst_pref = pref_soil_veg(jv,jst) DO ji=1,kjpindex corr_veg_soil(ji,jv,jst_pref) = zero corr_veg_soil_max(ji,jv,jst_pref) = zero !for veget distribution used in sechiba via humrel IF (soil_exist(ji,jst_pref).GT.min_sechiba) THEN corr_veg_soil(ji,jv,jst_pref)=MIN(veget_exist(ji,jv)/soiltype(ji,jst_pref),soil_exist(ji,jst_pref)) veget_exist(ji,jv)=MAX(veget_exist(ji,jv)-soil_exist(ji,jst_pref)*soiltype(ji,jst_pref),zero) soil_exist(ji,jst_pref)=MAX(soil_exist(ji,jst_pref)-corr_veg_soil(ji,jv,jst_pref),zero) ENDIF !same for max veget_max used in stomate via vegstress for slowproc IF (soil_exist_max(ji,jst_pref).GT.min_sechiba) THEN corr_veg_soil_max(ji,jv,jst_pref)= & & MIN(veget_exist_max(ji,jv)/soiltype(ji,jst_pref),soil_exist_max(ji,jst_pref)) veget_exist_max(ji,jv)=MAX(veget_exist_max(ji,jv)-soil_exist_max(ji,jst_pref)*soiltype(ji,jst_pref),zero) soil_exist_max(ji,jst_pref)=MAX(soil_exist_max(ji,jst_pref)-corr_veg_soil_max(ji,jv,jst_pref),zero) ENDIF ENDDO ENDDO ENDDO ENDIF ! ! update the corresponding masks ! ! mask_veget(:,:) = MIN( un, MAX(zero,veget(:,:))) ! mask_corr_veg_soil(:,:,:) = MIN( un, MAX(zero,corr_veg_soil(:,:,:))) mask_veget(:,:) = 0 mask_corr_veg_soil(:,:,:) = 0 DO ji = 1, kjpindex DO jv = 1, nvm IF(veget(ji,jv) .GT. min_sechiba) THEN mask_veget(ji,jv) = 1 ENDIF DO jst = 1, nstm IF(corr_veg_soil(ji,jv,jst) .GT. min_sechiba) THEN mask_corr_veg_soil(ji,jv,jst) = 1 ENDIF END DO END DO ! WRITE(numout,*) 'mask: soiltype,mask_soiltype',soiltype(ji,:),mask_soiltype(ji,:) END DO ! RETURN ! END SUBROUTINE hydrol_vegupd !! !! this routine computes soil processes with CWRR scheme !! SUBROUTINE hydrol_soil (kjpindex, dtradia, veget, veget_max, soiltype, transpir, vevapnu, evapot, & & evapot_penm, runoff, drainage, returnflow, irrigation, & & tot_melt, evap_bare_lim, shumdiag, litterhumdiag, humrel,vegstress, drysoil_frac) ! ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex ! input fields REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Map of vegetation types REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Map of max vegetation types REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: transpir !! transpiration REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapnu !! 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) :: evapot !! REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: evapot_penm !! ! modified fields REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: runoff !! complete runoff REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: drainage !! complete drainage REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_melt REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: evap_bare_lim !! 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,nvm), INTENT (inout) :: humrel !! 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) :: drysoil_frac !! Function of the litter humidity, that will be used to compute albedo ! ! local declaration ! INTEGER(i_std) :: ji, jv, jsl, jsl1, jst, ji_nsat !! indices INTEGER(i_std) :: m_sl0, m_sl1, isat !! mask values REAL(r_std) :: dztmp !! temporary depth REAL(r_std) :: temp !! temporary value for fluxes REAL(r_std) :: dpue !! temporary depth REAL(r_std), DIMENSION(kjpindex) :: tmcold, tmcint REAL(r_std), DIMENSION(kjpindex,nslm,nstm) :: moderwilt REAL(r_std), DIMENSION(kjpindex,nslm) :: mcint !! To save mc values for future use REAL(r_std), DIMENSION(kjpindex) :: correct_excess !! Corrects the flux at nslme layer in case of under residual moisture REAL(r_std), DIMENSION(kjpindex) :: mce !! Same use as mcint but jut for last efficient layer nslme REAL(r_std), DIMENSION(kjpindex) :: under_mcr !! Allows under residual soil moisture due to evap REAL(r_std), DIMENSION(kjpindex,nstm) :: v2 REAL(r_std), DIMENSION(kjpindex,nstm) :: evap_bare_lim_ns !! limitation of bare soi evaporation on each soil column (used to deconvoluate vevapnu) REAL(r_std) :: deltahum,diff REAL(r_std), DIMENSION(kjpindex) :: tsink REAL(r_std), DIMENSION(kjpindex) :: returnflow_soil REAL(r_std), DIMENSION(kjpindex) :: irrigation_soil REAL(r_std), DIMENSION(kjpindex,nstm) :: runoff_excess !! Runoff generated after soil saturation REAL(r_std) :: excess LOGICAL :: propagate !! if we propagate an excess ! ! returnflow_soil(:) = zero irrigation_soil(:) = zero qflux(:,:,:) = zero mask_return(:) = 0 index_nsat(:,:) = 0 index_sat(:,:) = 0 nslme(:,:) = nslm runoff_excess(:,:) = zero mce(:) = zero under_mcr(:) = zero correct_excess(:) = zero free_drain_coef(:,:) = zero n_sat(:) = 0 n_nsat(:) = 1 ! ! split 2d variables to 3d variables, per soil type ! CALL hydrol_split_soil (kjpindex, veget, soiltype, vevapnu, transpir, humrel, evap_bare_lim) ! ! for each soil type ! DO ji=1,kjpindex IF(vegtot(ji).GT.min_sechiba) THEN returnflow_soil(ji) = returnflow(ji)/vegtot(ji) irrigation_soil(ji) = irrigation(ji)/vegtot(ji) ENDIF DO jst=1, nstm !- A priori on considere qu'on est non sature si soiltype > 0 index_nsat(n_nsat(jst),jst)= ji*mask_soiltype(ji,jst) n_nsat(jst)=n_nsat(jst)+mask_soiltype(ji,jst) ENDDO IF(returnflow_soil(ji).GT.min_sechiba) THEN mask_return(ji) = 1 DO jst= 1,nstm isat=1 nslme(ji,jst)=nslm-isat ! DO jsl= nslm,3,-1 IF(mcs(jst)-mc(ji,jsl,jst).LT.min_sechiba) THEN nslme(ji,jst) = nslme(ji,jst) - isat ELSE isat = 0 ENDIF ENDDO ! We compute the indeces of the non-saturated points IF (nslme(ji,jst).LT.2) THEN nslme(ji,jst) = 0 ! En fait on est sature! n_nsat(jst)=n_nsat(jst)-1 n_sat(jst) = n_sat(jst)+1 index_sat(n_sat(jst),jst)=ji ENDIF ENDDO ENDIF ENDDO DO jst = 1,nstm n_nsat(jst) = n_nsat(jst)-1 ENDDO DO jst = 1,nstm ! !- We compute the sum of the sinks for future check-up DO ji=1,kjpindex tsink(ji) = SUM(rootsink(ji,:,jst))+MAX(ae_ns(ji,jst),zero)+subsinksoil(ji) ENDDO ! We save the Total moisture content tmcold(:) = tmc(:,jst) DO ji_nsat=1,n_nsat(jst) ji = index_nsat(ji_nsat,jst) !- the bare soil evaporation is substracted to the soil moisture profile: dpue = zz(nslme(ji,jst),jst) + dz(nslme(ji,jst)+1,jst) / deux DO jsl = 1, nslme(ji,jst) mc(ji,jsl,jst) = mc(ji,jsl,jst) & & - (MAX(ae_ns(ji,jst),zero) + subsinksoil(ji)) / dpue ENDDO !- we add the returnflow to the last efficient layer in the soil mc(ji,nslme(ji,jst),jst) = mc(ji,nslme(ji,jst),jst) + deux * returnflow_soil(ji) & & / (dz(nslme(ji,jst),jst) + dz(nslme(ji,jst)+1,jst)) ENDDO !!! SUBROUTINE hydrol_avoid_underres !-when mc(ji,1,jst)