Changeset 7709


Ignore:
Timestamp:
2022-07-20T11:30:43+02:00 (2 years ago)
Author:
josefine.ghattas
Message:

Integrated new irrigation scheme developed by Pedro Arboleda. See ticket #857
This corresponds to revsion 7708 of version pedro.arboleda/ORCHIDEE. Following differences were made but were not made on the pedro.arboleda/ORCHIDEE :

  • argumet place in call to routing_wrapper_intialize changed order
  • lines with only change in space were not taken
  • some indentation changed
  • set irrigation output as enalbled false if not do_irrigation
Location:
branches/ORCHIDEE_2_2/ORCHIDEE
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_parallel/xios_orchidee.f90

    r7576 r7709  
    451451          CALL xios_set_field_attr("netirrig",enabled=.FALSE.) 
    452452          CALL xios_set_field_attr("irrigmap",enabled=.FALSE.) 
     453          CALL xios_set_field_attr("irrig_deficit",enabled=.FALSE.) 
     454          CALL xios_set_field_attr("irrig_adduct",enabled=.FALSE.) 
     455          CALL xios_set_field_attr("irrig_gw_source",enabled=.FALSE.) 
     456          CALL xios_set_field_attr("irrig_fast_source",enabled=.FALSE.) 
     457          CALL xios_set_field_attr("irrig_str_source",enabled=.FALSE.) 
     458          CALL xios_set_field_attr("Count_failure_slow",enabled=.FALSE.) 
     459          CALL xios_set_field_attr("Count_failure_fast",enabled=.FALSE.) 
     460          CALL xios_set_field_attr("Count_failure_stre",enabled=.FALSE.) 
     461          CALL xios_set_field_attr("irrigmap_dyn",enabled=.FALSE.) 
    453462       END IF 
    454463 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes.f90

    r7547 r7709  
    879879 
    880880 
     881 
     882    !Crop irrigation 
     883    !Config Key   = IRRIG_DOSMAX 
     884    !Config Desc  = The maximum irrigation water injected per hour (kg.m^{-2}/hour) 
     885    !Config If    = DO_IRRIGATION AND old_irrig_scheme = FALSE 
     886    !Config Def   = 1. 
     887    !Config Help  = 
     888    !Config Units = [kg.m^{-2}/hour] 
     889    CALL getin_p('IRRIG_DOSMAX',irrig_dosmax) 
     890 
     891    !Config Key   = CUM_NROOT_THR 
     892    !Config Desc  = Cumulated nroot threshoold to define root zone, and calculate water deficit for irrigation 
     893    !Config If    = DO_IRRIGATION AND old_irrig_scheme = FALSE 
     894    !Config Def   = 0.90 
     895    !Config Help  = 
     896    !Config Units = [] 
     897    CALL getin_p('CUM_NROOT_THR',cum_nroot_thr) 
     898 
     899    !Crop irrigation 
     900    !Config Key   = IRRIGATED_SOILTILE 
     901    !Config Desc  = Do we introduce a new soil tile for irrigated croplands? 
     902    !Config If    = DO_IRRIGATION AND old_irrig_scheme = FALSE 
     903    !Config Def   = n 
     904    !Config Help  = 
     905    !Config Units = [FLAG] 
     906    CALL getin_p('IRRIGATED_SOILTILE',irrigated_soiltile) 
     907 
     908    ! 
     909    !Config Key   = OLD_IRRIG_SCHEME 
     910    !Config Desc  = Do we run with the old irrigation scheme? 
     911    !Config If    = DO_IRRIGATION 
     912    !Config Def   = n 
     913    !Config Help  = 
     914    !Config Units = [FLAG] 
     915    CALL getin_p('OLD_IRRIG_SCHEME',old_irrig_scheme) 
     916 
     917    !   To run with new irrigation scheme 
     918    !Config Key   = IRRIG_ST 
     919    !Config Desc  = Which is the soil tile with irrigation flux 
     920    !Config If    = DO_IRRIGATION and old_irrig_scheme = FALSE and IRRIGATED_SOILTILE = FALSE 
     921    !Config Def   = 3 
     922    !Config Help  = 
     923    !Config Units = [] 
     924    CALL getin_p('IRRIG_ST',irrig_st) 
     925 
     926    !  To run with NEW irrigation scheme 
     927    !Config Key   = AVAIL_RESERVE 
     928    !Config Desc  = Max. fraction of routing reservoir volume that can be used for irrigation 
     929    !Config If    = DO_IRRIGATION and old_irrig_scheme = FALSE 
     930    !Config Def   = 0.9,0.0,0.9 
     931    !Config Help  = 
     932    !Config Units = [] 
     933    CALL getin_p('AVAIL_RESERVE',avail_reserve) 
     934 
     935    !  To run with NEW irrigation scheme 
     936    !Config Key   = BETA_IRRIG 
     937    !Config Desc  = Threshold multiplier of Target SM to calculate root deficit 
     938    !Config If    = DO_IRRIGATION and old_irrig_scheme = FALSE 
     939    !Config Def   = 1. 
     940    !Config Help  = 
     941    !Config Units = [] 
     942    CALL getin_p('BETA_IRRIG',beta_irrig) 
     943    ! 
     944    !   To run with new irrigation scheme 
     945    !Config Key   = LAI_IRRIG_MIN 
     946    !Config Desc  = Min. LAI (Leaf Area Index) to trigger irrigation 
     947    !Config If    = DO_IRRIGATION and old_irrig_scheme = FALSE and IRRIGATED_SOILTILE = FALSE 
     948    !Config Def   = 0.1 
     949    !Config Help  = 
     950    !Config Units = [m2/m2] 
     951    CALL getin_p('LAI_IRRIG_MIN',lai_irrig_min) 
     952 
     953    !   To run with new irrigation scheme 
     954    !Config Key   = irrig_map_dynamic 
     955    !Config Desc  = Do we use a dynamic irrig map? 
     956    !Config If    = DO_IRRIGATION = T 
     957    !Config Def   = n 
     958    !Config Help  = 
     959    !Config Units = [FLAG] 
     960    CALL getin_p('IRRIG_MAP_DYNAMIC_FLAG',irrig_map_dynamic_flag) 
     961 
     962    !   To run with new irrigation scheme 
     963    !Config Key   = select_source_irrig 
     964    !Config Desc  = Do we use the new priorization scheme, based on maps of equipped area with surface water? 
     965    !Config If    = DO_IRRIGATION = T 
     966    !Config Def   = n 
     967    !Config Help  = 
     968    !Config Units = [FLAG] 
     969    CALL getin_p('SELECT_SOURCE_IRRIG',select_source_irrig) 
     970 
     971    !   Do we reinfiltrate all/part of runoff in crop soil tile? 
     972    !Config Key   = Reinfiltr_IrrigField 
     973    !Config Desc  = DO_IRRIGATION 
     974    !Config If    = DO_IRRIGATION = T 
     975    !Config Def   = n 
     976    !Config Help  = 
     977    !Config Units = [FLAG] 
     978    CALL getin_p('REINFILTR_IRRIGFIELD',Reinfiltr_IrrigField) 
     979 
     980    !   Externalized reinf_slope for reinfiltration in irrigated cropland 
     981    !Config Key   = reinf_slope_cropParam 
     982    !Config Desc  = DO_IRRIGATION 
     983    !Config If    = DO_IRRIGATION = T, REINFILTR_IRRIGFIELD = T 
     984    !Config Def   = 0.8 
     985    !Config Help  = 
     986    !Config Units = [FRACTION 0-1] 
     987    CALL getin_p('REINF_SLOPE_CROPPARAM',reinf_slope_cropParam) 
     988 
     989    !   Externalized for available volume to adduction 
     990    !Config Key   = a_stream_adduction   
     991    !Config Desc  = DO_IRRIGATION 
     992    !Config If    = DO_IRRIGATION = T 
     993    !Config Def   = zero 
     994    !Config Help  = 
     995    !Config Units = [FRACTION 0-1] 
     996    CALL getin_p('A_STREAM_ADDUCTION',a_stream_adduction) 
     997 
     998 
    881999    !Surface resistance 
    8821000    ! 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_var.f90

    r7547 r7709  
    497497!$OMP THREADPRIVATE(methanol_activity) 
    498498 
     499 
     500 ! 
     501 ! Parameters for irrigation scheme 
     502 ! 
     503  REAL(r_std), SAVE :: irrig_dosmax = 1.              !! The maximum irrigation water injected per hour (kg.m^{-2}/hour) 
     504!$OMP THREADPRIVATE(irrig_dosmax) 
     505  REAL(r_std), SAVE :: cum_nroot_thr = 0.90           !! Cumulated nroot threshoold to define root zone, and calculate water deficit for irrigation (-) 
     506!$OMP THREADPRIVATE(cum_nroot_thr) 
     507  LOGICAL, SAVE :: irrigated_soiltile = .FALSE.       !! Do we introduce a new soil tile for irrigated croplands? (true/false) 
     508!$OMP THREADPRIVATE(irrigated_soiltile) 
     509  LOGICAL, SAVE :: old_irrig_scheme = .FALSE.         !! Do we run with the old irrigation scheme? (true/false)  , add to compatiblity 
     510!$OMP THREADPRIVATE(old_irrig_scheme) 
     511  INTEGER, SAVE :: irrig_st = 3                       !! Which is the soil tile with irrigation flux 
     512!$OMP THREADPRIVATE(irrig_st) 
     513  REAL(r_std), SAVE, DIMENSION(3) :: avail_reserve = (/0.9,0.0,0.9/)     !! Available water from routing reservoirs, to withdraw for irrigation 
     514                                                      !! IMPORTANT: As the routing model uses 3 reservoirs, dimension is set to 3 
     515                                                      !! IMPORTANT: Order of available water must be in this order: streamflow, fast, and slow reservoir 
     516!$OMP THREADPRIVATE(avail_reserve) 
     517  REAL(r_std), SAVE :: beta_irrig = 1.                !! Threshold multiplier of Target SM to calculate root deficit(unitless) 
     518!$OMP THREADPRIVATE(beta_irrig) 
     519  REAL(r_std), SAVE :: lai_irrig_min = 0.1            !! Minimum LAI to trigger irrigation (kg.m^{-2}/hour) 
     520!$OMP THREADPRIVATE(lai_irrig_min) 
     521  LOGICAL, SAVE :: irrig_map_dynamic_flag = .FALSE.   !! Do we use a dynamic irrig map? 
     522!$OMP THREADPRIVATE(irrig_map_dynamic_flag) 
     523  LOGICAL, SAVE :: select_source_irrig = .FALSE.      !! Do we use the new priorization scheme, based on maps of equipped area with surface water? 
     524!$OMP THREADPRIVATE(select_source_irrig) 
     525  LOGICAL, SAVE :: Reinfiltr_IrrigField = .FALSE.     !! Do we reinfiltrate all runoff from crop soil tile?O 
     526!$OMP THREADPRIVATE(Reinfiltr_IrrigField) 
     527  REAL, SAVE :: reinf_slope_cropParam = 0.8           !! Externalized for irrigated cropland, when Reinfiltr_IrrigField=.TRUE. 
     528                                                      !! Max value of reinf_slope in irrig_st   
     529!$OMP THREADPRIVATE(reinf_slope_cropParam) 
     530  REAL, SAVE :: a_stream_adduction = zero             !! Externalized for available volume to adduction 
     531!$OMP THREADPRIVATE(a_stream_adduction) 
     532 
     533 
     534 
    499535  ! 
    500536  ! condveg.f90 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/diffuco.f90

    r7326 r7709  
    332332    REAL(r_std),DIMENSION (kjpindex,nvm)     :: cim                !! Intercellular CO2 over nlai  
    333333 
     334    !! Based on ZunYin for correct transport 
     335    REAL(r_std),DIMENSION (kjpindex,nvm)     :: humrel_tmp         !! Soil moisture stress (within range 0 to 1) 
     336    REAL(r_std),DIMENSION (kjpindex,nvm)     :: vbeta3_tmp         !! Beta for transpiration (-) 
     337    REAL(r_std),DIMENSION (kjpindex,nvm)     :: rveget_tmp         !! Stomatal resistance for the whole canopy (s m^{-1}) 
     338    REAL(r_std),DIMENSION (kjpindex,nvm)     :: rstruct_tmp        !! Structural resistance for the vegetation 
     339    REAL(r_std),DIMENSION (kjpindex,nvm)     :: cimean_tmp         !! Mean leaf Ci (ppm) 
     340    REAL(r_std),DIMENSION (kjpindex,nvm)     :: gsmean_tmp         !! Mean stomatal conductance to CO2 (umol m-2 s-1) 
     341    REAL(r_std),DIMENSION (kjpindex,nvm)     :: gpp_tmp            !! Assimilation ((gC m^{-2} s^{-1}), total area) 
     342    REAL(r_std),DIMENSION (kjpindex,nvm)     :: cim_tmp            !! Intercellular CO2 over nlai 
     343 
    334344!_ ================================================================================================================================ 
    335345     
    … …  
    368378 
    369379  !! 5. Calculate the beta coefficient for transpiration 
     380    !! Dummy variables to compute vbeta3pot (humrel_temp = 1, no water stress) 
     381    humrel_tmp(:,:) = un 
     382    vbeta3_tmp = vbeta3 
     383    rveget_tmp = rveget 
     384    rstruct_tmp = rstruct 
     385    cimean_tmp = cimean 
     386    gsmean_tmp = gsmean 
     387    gpp_tmp = gpp 
     388    cim_tmp = cim 
    370389 
    371390    CALL diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel, & 
    … …  
    374393         rveget, rstruct, cimean, gsmean, gpp, & 
    375394         vbeta23, hist_id, indexveg, indexlai, index, kjit, cim) 
     395 
     396   !! New resistance calc. with temporary (dummy) variables, just for vbeta3pot 
     397    CALL diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel_tmp, & 
     398         assim_param, ccanopy, & 
     399         veget, veget_max, lai, qsintveg, qsintmax, vbeta3_tmp, vbeta3pot, & 
     400         rveget_tmp, rstruct_tmp, cimean_tmp, gsmean_tmp, gpp_tmp, & 
     401         vbeta23, hist_id, indexveg, indexlai, index, kjit, cim_tmp) 
    376402 
    377403    ! 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90

    r7476 r7709  
    2424!!                    We can also impose kfact_root=1 in all soil layers to cancel the effect of 
    2525!!                    roots on ks profile (keyword KFACT_ROOT_CONST). 
     26!!                    July 2022: New irrigation scheme. Here new irrigation demand based in 
     27!!                    soil moisture deficit, and irrigation application. 
    2628!!                  
    2729!! REFERENCE(S) : 
    … …  
    366368                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex 
    367369!$OMP THREADPRIVATE(mc) 
     370  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: root_mc_fc         !! Max Field Capacity moisture content, for layers with roots, in soil tile of irrig_st 
     371                                                                         !!  @tex $(kg m^{-2})$ @endtex 
     372!$OMP THREADPRIVATE(root_mc_fc) 
     373  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: nslm_root          !! max. layers defining the root zone 
     374                                                                         !!  @tex $(layer)$ @endtex 
     375!$OMP THREADPRIVATE(nslm_root) 
    368376 
    369377   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc_read_prev       !! Soil moisture from file at previous timestep in the file 
    … …  
    565573       & humrel, vegstress, drysoil_frac, evapot, evapot_penm, evap_bare_lim, evap_bare_lim_ns, & 
    566574       & flood_frac, flood_res, & 
    567        & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope, rest_id, hist_id, hist2_id,& 
     575       & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope_soil, rest_id, hist_id, hist2_id,& 
    568576       & contfrac, stempdiag, & 
    569577       & temp_air, pb, u, v, tq_cdrag, swnet, pgflux, & 
    … …  
    571579       & grndflux,gtemp,tot_bare_soil, & 
    572580       & lambda_snow,cgrnd_snow,dgrnd_snow,frac_snow_veg,temp_sol_add, & 
    573        & mc_layh, mcl_layh, soilmoist_out ) 
     581       & mc_layh, mcl_layh, soilmoist_out, root_deficit ) 
    574582 
    575583    !! 0. Variable and parameter declaration 
    … …  
    607615    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: qsintmax         !! Maximum water on vegetation for interception 
    608616    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: transpir         !! Transpiration 
    609     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: reinf_slope      !! Slope coef 
     617    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: reinf_slope_soil !! Slope coef per soil tile 
    610618 
    611619    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1}) 
    … …  
    647655    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: tot_melt         !! Total melt     
    648656    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: floodout         !! Flux out of floodplains 
    649      
     657    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: root_deficit     !! water deficit to reach field capacity of soil 
     658 
    650659    !! 0.3 Modified variables 
    651660 
    … …  
    794803    !! 3.5 computes soil hydrology ==>hydrol_soil 
    795804 
    796     CALL hydrol_soil(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, reinf_slope,  & 
     805    CALL hydrol_soil(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, reinf_slope_soil,  & 
    797806         transpir, vevapnu, evapot, evapot_penm, runoff, drainage, &  
    798807         returnflow, reinfiltration, irrigation, & 
    … …  
    800809         k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,& 
    801810         stempdiag,snow,snowdz, tot_bare_soil,  u, v, tq_cdrag, & 
    802          mc_layh, mcl_layh) 
     811         mc_layh, mcl_layh, root_deficit, veget) 
    803812 
    804813    ! The update fluxes come from hydrol_vegupd 
    805814    drainage(:) =  drainage(:) +  drain_upd(:) 
    806815    runoff(:) =  runoff(:) +  runoff_upd(:) 
     816 
     817   ! Output irrigation related variables 
     818    CALL xios_orchidee_send_field("root_deficit", root_deficit) 
     819    CALL xios_orchidee_send_field("root_mc_fc", root_mc_fc) 
    807820 
    808821 
    … …  
    16911704    ALLOCATE (mc(kjpindex,nslm,nstm),stat=ier) 
    16921705    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc','','') 
     1706 
     1707    ALLOCATE (root_mc_fc(kjpindex),stat=ier) ! 
     1708    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable root_mc_fc','','') ! 
     1709 
     1710    ALLOCATE (nslm_root(kjpindex),stat=ier) ! 
     1711    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nslm_root','','') ! 
    16931712 
    16941713 
    … …  
    23652384    IF (ALLOCATED  (vegetmax_soil)) DEALLOCATE (vegetmax_soil) 
    23662385    IF (ALLOCATED  (mc)) DEALLOCATE (mc) 
     2386    IF (ALLOCATED  (root_mc_fc)) DEALLOCATE (root_mc_fc) 
     2387    IF (ALLOCATED  (nslm_root)) DEALLOCATE (nslm_root) 
    23672388    IF (ALLOCATED  (soilmoist)) DEALLOCATE (soilmoist) 
    23682389    IF (ALLOCATED  (soilmoist_liquid)) DEALLOCATE (soilmoist_liquid) 
    … …  
    27282749    INTEGER(i_std)                                      :: iiref           !! To identify the mc_lins where k_lin and d_lin 
    27292750                                                                           !! need special treatment 
     2751    REAL(r_std), DIMENSION (kjpindex,nvm)               :: cum_nroot       !! local variable to acummulate nroot 
     2752    INTEGER(i_std), DIMENSION (kjpindex,nvm)            :: nslm_pft_root   !! Deeper layer with cum_nroot > cum_nroot_thr per PFT, 
     2753    REAL(r_std), DIMENSION (kjpindex,nslm)              :: smf             !! Soil moisture of each layer at field capacity 
     2754                                                                           !!  @tex $(kg m^{-2})$ @endtex 
    27302755 
    27312756!_ ================================================================================================================================ 
    … …  
    28802905    !! 2 Compute the root density profile if not ok_dynroot 
    28812906    !!   For the case with ok_dynroot, the calculations are done at each time step in hydrol_soil 
     2907    cum_nroot(:,:) = zero 
     2908    nslm_pft_root(:,:) = nslm 
     2909    nslm_root(:) =  nslm 
    28822910    IF (.NOT. ok_dynroot) THEN 
    2883        DO ji=1, kjpindex 
    2884           !- 
    2885           !! The three following equations concerning nroot computation are derived from the integrals  
    2886           !! of equations C9 to C11 of De Rosnay's (1999) PhD thesis (page 158). 
    2887           !! The occasional absence of minus sign before humcste parameter is correct. 
     2911       !! Calculation of nroot and of new nslm_root for irrigation  
     2912       !! The three following equations concerning nroot computation are derived from the integrals 
     2913       !! of equations C9 to C11 of De Rosnay's (1999) PhD thesis (page 158). 
     2914       !! The occasional absence of minus sign before humcste parameter is correct. 
     2915       ! First layer 
     2916       nroot(:,:,1) = zero 
     2917       !From 2 to nslm-1 layers 
     2918       DO jsl = 2, nslm-1 
    28882919          DO jv = 1,nvm 
    2889              DO jsl = 2, nslm-1 
     2920             DO ji=1, kjpindex 
     2921                 
     2922                 
    28902923                nroot(ji,jv,jsl) = (EXP(-humcste(jv)*zz(jsl)/mille)) * & 
    28912924                     & (EXP(humcste(jv)*dz(jsl)/mille/deux) - & 
    … …  
    28932926                     & (EXP(-humcste(jv)*dz(2)/mille/deux) & 
    28942927                     & -EXP(-humcste(jv)*zz(nslm)/mille)) 
     2928                !We acum. nroot, and change nslm_pft_root if necessary 
     2929                cum_nroot(ji,jv) = cum_nroot(ji,jv) + nroot(ji,jv,jsl) ! 
     2930                 
     2931                IF(cum_nroot(ji,jv) > cum_nroot_thr .AND. nslm_pft_root(ji,jv) .GE. jsl ) THEN 
     2932                   nslm_pft_root(ji,jv) =  jsl 
     2933                ENDIF 
    28952934             ENDDO 
    2896              nroot(ji,jv,1) = zero 
    2897  
     2935              
     2936          ENDDO 
     2937           
     2938       ENDDO 
     2939       !Last layer 
     2940       !No need to put and IF here, if it is the case, nslm_pft_root is already 
     2941       ! equal to nslm, 
     2942       DO jv = 1,nvm 
     2943          DO ji=1, kjpindex 
    28982944             nroot(ji,jv,nslm) = (EXP(humcste(jv)*dz(nslm)/mille/deux) -un) * & 
    28992945                  & EXP(-humcste(jv)*zz(nslm)/mille) / & 
    … …  
    29022948          ENDDO 
    29032949       ENDDO 
     2950       !New loop to compute min. soil layer per cell, using just PFT that 
     2951       !are not natural and exist inside the cell 
     2952       DO jv=1, nvm 
     2953          IF(.NOT. natural(jv)) THEN 
     2954             DO ji=1,kjpindex 
     2955                IF(veget_max(ji, jv) > min_sechiba) THEN 
     2956                   nslm_root(ji) = MIN(nslm_root(ji), nslm_pft_root(ji,jv)) 
     2957                ENDIF 
     2958             ENDDO 
     2959          ENDIF 
     2960       ENDDO 
     2961        
    29042962    END IF 
    29052963 
    2906   
    2907  
     2964    ! Calculates field capacity soil moisture per soil layers 
     2965    ! then calculate field capacity soil moisture over root zone 
     2966    smf(:,:) = zero 
     2967    root_mc_fc(:) = zero 
     2968    smf(:,1) = dz(2) * (quatre*mcfc(:))/huit 
     2969     
     2970    DO jsl = 2,nslm-1 
     2971       smf(:,jsl) = dz(jsl) * ( quatre*mcfc(:) )/huit & 
     2972            + dz(jsl+1) * ( quatre*mcfc(:) )/huit 
     2973    ENDDO 
     2974     
     2975    smf(:,nslm) = dz(nslm) * (quatre*mcfc(:))/huit 
     2976    DO ji = 1,kjpindex 
     2977       root_mc_fc(ji) = SUM(smf(ji,1:nslm_root(ji) )) 
     2978    ENDDO 
     2979     
    29082980    !- 
    29092981    !! 3 Compute the profile for a and n 
    … …  
    35923664!_ hydrol_soil 
    35933665  SUBROUTINE hydrol_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, & 
    3594        kjpindex, veget_max, soiltile, njsc, reinf_slope, & 
     3666       kjpindex, veget_max, soiltile, njsc, reinf_slope_soil, & 
    35953667       & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 
    35963668       & returnflow, reinfiltration, irrigation, & 
    … …  
    35983670       & k_litt, litterhumdiag, humrel,vegstress, drysoil_frac, & 
    35993671       & stempdiag,snow, & 
    3600        & snowdz, tot_bare_soil, u, v, tq_cdrag, mc_layh, mcl_layh) 
     3672       & snowdz, tot_bare_soil, u, v, tq_cdrag, mc_layh, mcl_layh, root_deficit, veget) 
    36013673    !  
    36023674    ! interface description 
    … …  
    36083680    INTEGER(i_std), INTENT(in)                               :: kjpindex  
    36093681    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in)       :: veget_max        !! Map of max vegetation types [-] 
     3682    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: veget            !! Fraction of vegetation type 
    36103683    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: njsc             !! Index of the dominant soil textural class  
    36113684                                                                                 !!   in the grid cell (1-nscm, unitless) 
    … …  
    36223695    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: transpir         !! Transpiration   
    36233696                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex 
    3624     REAL(r_std), DIMENSION (kjpindex), INTENT (in)           :: reinf_slope      !! Fraction of surface runoff that reinfiltrates 
     3697    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in)      :: reinf_slope_soil !! Fraction of surface runoff that reinfiltrates per soil tile 
    36253698                                                                                 !!  (unitless, [0-1]) 
    36263699    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: returnflow       !! Water returning to the soil from the bottom 
    … …  
    36723745    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: mcl_layh         !! Volumetric liquid water content for each soil layer  
    36733746                                                                                 !! averaged over the mesh (for thermosoil) 
    3674                                                                                  !!  @tex $(m^{3} m^{-3})$ @endtex  
     3747                                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex      
     3748    REAL(r_std),DIMENSION (kjpindex), INTENT(out)            :: root_deficit     !! water deficit to reach SM target of soil column, for irrigation demand 
     3749     
    36753750    !! 0.3 Modified variables 
    36763751 
    … …  
    36983773    REAL(r_std), DIMENSION(kjpindex)               :: reinfiltration_soil        !! Water from the routing back to the top of the  
    36993774                                                                                 !! soil @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex 
    3700     REAL(r_std), DIMENSION(kjpindex)               :: irrigation_soil            !! Water from irrigation returning to soil moisture  
     3775    REAL(r_std), DIMENSION(kjpindex,nstm)          :: irrigation_soil            !! Water from irrigation returning to soil moisture per soil tile 
    37013776                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex 
    37023777    REAL(r_std), DIMENSION(kjpindex)               :: flux_infilt                !! Water to infiltrate 
    … …  
    37663841    REAL(r_std), DIMENSION (kjpindex)              :: check_over                  !! Water conservation diagnostic at routine scale  
    37673842    REAL(r_std), DIMENSION (kjpindex)              :: check_under                 !! Water conservation diagnostic at routine scale 
    3768  
     3843    ! For irrigation triggering 
     3844    INTEGER(i_std), DIMENSION (kjpindex)           :: lai_irrig_trig              !! Number of PFT per cell with LAI> LAI_IRRIG_MIN - 
    37693845    ! Diagnostic of the vertical soil water fluxes   
    37703846    REAL(r_std), DIMENSION (kjpindex,nslm)         :: qflux                       !! Local upward flux into soil layer  
    … …  
    37963872    returnflow_soil(:) = zero 
    37973873    reinfiltration_soil(:) = zero 
    3798     irrigation_soil(:) = zero 
     3874    irrigation_soil(:,:) = zero 
    37993875    qflux_ns(:,:,:) = zero 
    38003876    mc_layh(:,:) = zero ! for thermosoil 
    … …  
    38383914    ! AD16*** The transformation in 0.2 and 0.3 is likely to induce conservation problems 
    38393915    !         when tot_frac_nobio NE 0, since sum(soiltile) NE vegtot in this case 
    3840      
    3841     DO ji=1,kjpindex 
    3842        IF(vegtot(ji).GT.min_sechiba) THEN 
    3843           ! returnflow_soil is assumed to enter from the bottom, but it is not possible with CWRR 
    3844           returnflow_soil(ji) = zero 
    3845           reinfiltration_soil(ji) = (returnflow(ji) + reinfiltration(ji))/vegtot(ji) 
    3846           irrigation_soil(ji) = irrigation(ji)/vegtot(ji) 
    3847        ELSE 
    3848           returnflow_soil(ji) = zero 
    3849           reinfiltration_soil(ji) = zero 
    3850           irrigation_soil(ji) = zero 
     3916    IF (.NOT. old_irrig_scheme) THEN ! 
     3917       IF (.NOT. irrigated_soiltile) THEN 
     3918          DO ji=1,kjpindex 
     3919             IF(vegtot(ji).GT.min_sechiba ) THEN 
     3920                returnflow_soil(ji) = zero 
     3921                reinfiltration_soil(ji) = (returnflow(ji) + reinfiltration(ji))/vegtot(ji) 
     3922                IF(soiltile(ji, irrig_st).GT.min_sechiba) THEN 
     3923                   !irrigation_soil(ji, 1:2) = 0, if irrig_st = 3. Not put because Values 
     3924                   !are already zero, due to initialization 
     3925                   irrigation_soil(ji, irrig_st) = irrigation(ji) / (soiltile(ji, irrig_st) * vegtot(ji) ) 
     3926                   !Irrigation is kg/m2 of grid cell. Here, all that water is put on 
     3927                   !irrig_st (irrigated soil tile), by default = 3, for the others 
     3928                   !soil tiles irrigation = zero 
     3929                ENDIF 
     3930             ENDIF 
     3931          ENDDO 
    38513932       ENDIF 
    3852     ENDDO        
     3933    ELSE  ! 
     3934       DO ji=1,kjpindex 
     3935          IF(vegtot(ji).GT.min_sechiba) THEN 
     3936             ! returnflow_soil is assumed to enter from the bottom, but it is not possible with CWRR 
     3937             returnflow_soil(ji) = zero 
     3938             reinfiltration_soil(ji) = (returnflow(ji) + reinfiltration(ji))/vegtot(ji) 
     3939             irrigation_soil(ji,:) = irrigation(ji)/vegtot(ji) 
     3940             ! irrigation_soil(ji) = irrigation(ji)/vegtot(ji) 
     3941             ! Computed for all the grid cell. New way is equivalent, and coherent 
     3942             ! with irrigation_soil new dimensions (cells, soil tiles) 
     3943             ! Irrigation is kg/m2 of grid cell. For the old irrig. scheme, 
     3944             ! irrigation soil is the same for every soil tile 
     3945             ! Next lines are in tag 2.0, deleted because values are already init to zero 
     3946             ! ELSE 
     3947             ! returnflow_soil(ji) = zero 
     3948             ! reinfiltration_soil(ji) = zero 
     3949             ! irrigation_soil(ji) = zero 
     3950             ! ENDIF 
     3951          ENDIF 
     3952       ENDDO 
     3953    ENDIF 
    38533954     
    38543955    !! -- START MAIN LOOP (prognostic loop to update mc and mcl) OVER SOILTILES 
    … …  
    39014002       ! Here, temp is used as a temporary variable to store the min of water to infiltrate vs evaporate 
    39024003       DO ji = 1, kjpindex 
    3903           temp(ji) = MIN(water2infilt(ji,jst) + irrigation_soil(ji) + reinfiltration_soil(ji) & 
     4004          temp(ji) = MIN(water2infilt(ji,jst) + irrigation_soil(ji,jst) + reinfiltration_soil(ji) & 
    39044005                         - MIN(ae_ns(ji,jst),zero) - MIN(subsinksoil(ji),zero) + precisol_ns(ji,jst), & 
    39054006                           MAX(ae_ns(ji,jst),zero) + MAX(subsinksoil(ji),zero) ) 
    … …  
    39124013       !   - eventually, water2infilt holds all fluxes to the soil surface except precisol (reduced by water2extract) 
    39134014       DO ji = 1, kjpindex 
    3914           water2infilt(ji,jst) = water2infilt(ji,jst) + irrigation_soil(ji) + reinfiltration_soil(ji) & 
     4015         !Note that in tag 2.0, irrigation_soil(ji), changed to be coherent with new variable dimension 
     4016          water2infilt(ji,jst) = water2infilt(ji,jst) + irrigation_soil(ji, jst) + reinfiltration_soil(ji) & 
    39154017                - MIN(ae_ns(ji,jst),zero) - MIN(subsinksoil(ji),zero) + precisol_ns(ji,jst) & 
    39164018                - temp(ji)  
    … …  
    39554057       IF ( .NOT. doponds ) THEN ! this is the general case... 
    39564058          DO ji = 1, kjpindex 
    3957              water2infilt(ji,jst) = reinf_slope(ji) * ru_ns(ji,jst) 
     4059             water2infilt(ji,jst) = reinf_slope_soil(ji,jst) * ru_ns(ji,jst) 
    39584060          ENDDO 
    39594061       ELSE 
    … …  
    43414443          tmcs_litter(:) = tmcs_litter(:) + sms(:,jsl) 
    43424444       END DO 
    4343               
     4445 
     4446       ! Here we compute root zone deficit, to have an estimate of water demand in irrigated soil column (i.e. crop and grass) 
     4447       IF(jst .EQ. irrig_st ) THEN 
     4448           !It computes water deficit only on the root zone, and only on the layers where 
     4449           !there is actually a deficit. If there is not deficit, it does not take into account that layer 
     4450            DO ji = 1,kjpindex 
     4451 
     4452                root_deficit(ji) = SUM( MAX(zero, beta_irrig*smf(ji,1:nslm_root(ji) ) & 
     4453                - sm(ji,1:nslm_root(ji)  ) )) - water2infilt(ji,jst) 
     4454 
     4455                root_deficit(ji) = MAX( root_deficit(ji) ,  zero) 
     4456 
     4457            ENDDO 
     4458            !It COUNTS the number of pft with LAI > lai_irrig_min, inside the soiltile 
     4459            !It compares veget, but it is the same as they are related by a function 
     4460            lai_irrig_trig(:) = 0 
     4461 
     4462            DO jv = 1, nvm 
     4463              IF( .NOT. natural(jv) ) THEN 
     4464                DO ji = 1,kjpindex 
     4465 
     4466                  IF(veget(ji, jv) > veget_max(ji, jv) * ( un - & 
     4467                  exp( - lai_irrig_min * ext_coeff_vegetfrac(jv) ) ) ) THEN 
     4468 
     4469                    lai_irrig_trig(ji) = lai_irrig_trig(ji) + 1 
     4470 
     4471                  ENDIF 
     4472 
     4473                ENDDO 
     4474              ENDIF 
     4475 
     4476            ENDDO 
     4477            !If any of the PFT inside the soil tile have LAI >  lai_irrig_min (I.E. lai_irrig_trig(ji) = 0 ) 
     4478            !The root deficit is set to zero, and irrigation is not triggered 
     4479            DO ji = 1,kjpindex 
     4480 
     4481              IF( lai_irrig_trig(ji) < 1 ) THEN 
     4482                root_deficit(ji) = zero 
     4483              ENDIF 
     4484 
     4485            ENDDO 
     4486       ENDIF 
     4487 
    43444488       ! Soil wetness profiles (W-Ww)/(Ws-Ww) 
    43454489       ! soil_wet_ns is the ratio of available soil moisture to max available soil moisture 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90

    r7576 r7709  
    1212!!\n DESCRIPTION: None 
    1313!! 
    14 !! RECENT CHANGE(S): None 
    15 !! 
     14!! RECENT CHANGE(S): July 2022: New irrigation scheme. Here new irrigation offer with more information from maps, and 
     15!!                    calculation of water withdrawal. 
    1616!! REFERENCE(S) : 
    1717!! 
    … …  
    226226  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: reinfiltration_mean          !! Mean water flow which returns to the grid box (kg/m^2/dt) 
    227227!$OMP THREADPRIVATE(reinfiltration_mean) 
     228  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrigdeficit_mean              !! Mean irrigation deficit. 
     229                                                                                           !! This is between irrigation and water requirement (kg/m^2/dt) 
     230!$OMP THREADPRIVATE(irrigdeficit_mean) 
    228231  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrigation_mean              !! Mean irrigation flux. 
    229232                                                                                             !! This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt) 
    230233!$OMP THREADPRIVATE(irrigation_mean) 
     234  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrig_gw_source_mean       !! Mean groundwater irrigation flux. 
     235                                                                                         !! This is the water taken from the GW reservoir only (kg/m^2/dt) 
     236!$OMP THREADPRIVATE(irrig_gw_source_mean) 
     237  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrig_fast_source_mean       !! Mean irrigation flux from fast. 
     238                                                                                       !! This is the water taken from the fast reservoir only (kg/m^2/dt) 
     239!$OMP THREADPRIVATE(irrig_fast_source_mean) 
     240  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrig_str_source_mean       !! Mean streamflow irrigation flux. 
     241                                                                                       !! This is the water taken from the streamflow reservoir only (kg/m^2/dt) 
     242!$OMP THREADPRIVATE(irrig_str_source_mean) 
     243  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrigadduct_mean       !! Irrigation that comes from adduction. It includes water from basins inside the grid cell 
     244                                                                                        ! and water from nearby grid cells 
     245!$OMP THREADPRIVATE(irrigadduct_mean) 
    231246  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: riverflow_mean               !! Mean Outflow of the major rivers. 
    232247                                                                                             !! The flux will be located on the continental grid but this should be a coastal point (kg/dt) 
    … …  
    291306                                neighbours,  resolution,     contfrac,   stempdiag, & 
    292307                                returnflow,  reinfiltration, irrigation, riverflow, & 
    293                                 coastalflow, flood_frac,     flood_res ) 
    294         
     308                                coastalflow, flood_frac,     flood_res, soiltile, irrig_frac, & 
     309                                veget_max, irrigated_next) ! 
     310 
    295311    IMPLICIT NONE 
    296312     
    … …  
    309325    REAL(r_std), INTENT(in)        :: contfrac(nbpt)       !! Fraction of land in each grid box (unitless;0-1) 
    310326    REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile 
     327    REAL(r_std), INTENT(in)        :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless) 
     328    REAL(r_std), INTENT(in)        :: veget_max(nbpt,nvm)  !! Maximal fraction of vegetation (unitless;0-1) ! 
     329    REAL(r_std), INTENT(in)        :: irrigated_next (nbpt)  !! Dynamic irrig. area, calculated in slowproc and passed to routing! 
     330    REAL(r_std), INTENT(in)        :: irrig_frac(nbpt)      !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE. 
    311331 
    312332    !! 0.2 Output variables 
    … …  
    323343    !! 0.3 Local variables 
    324344    REAL(r_std), DIMENSION(nbp_glo):: mask_coast_glo       !! Mask with coastal gridcells on global grid (1/0) 
    325     LOGICAL                        :: init_irrig           !! Logical to initialize the irrigation (true/false) 
     345    !LOGICAL                        :: init_irrig          !! Logical to initialize the irrigation (true/false) 
    326346    LOGICAL                        :: init_flood           !! Logical to initialize the floodplains (true/false) 
    327347    LOGICAL                        :: init_swamp           !! Logical to initialize the swamps (true/false) 
    … …  
    385405    !! Initialisation of flags for irrigated land, flood plains and swamps 
    386406    ! 
    387     init_irrig = .FALSE. 
    388407    IF ( do_irrigation ) THEN  
    389        IF (COUNT(irrigated .GE. undef_sechiba-1) > 0) init_irrig = .TRUE. 
     408       irrigated(:) = irrigated_next(:) 
    390409    END IF 
    391410     
    … …  
    399418       IF (COUNT(swamp .GE. undef_sechiba-1) > 0 ) init_swamp = .TRUE. 
    400419    END IF 
    401         
     420 
    402421    !! If we have irrigated land, flood plains or swamps then we need to interpolate the 0.5 degree 
    403422    !! base data set to the resolution of the model. 
    404      
    405     IF ( init_irrig .OR. init_flood .OR. init_swamp ) THEN 
    406        CALL routing_irrigmap(nbpt, index, lalo, neighbours, resolution, & 
    407             contfrac, init_irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id) 
     423    !IF ( init_irrig .OR. init_flood .OR. init_swamp ) THEN 
     424    IF ( init_flood .OR. init_swamp ) THEN 
     425       CALL routing_floodmap(nbpt, index, lalo, neighbours, resolution, & 
     426            contfrac, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id) 
    408427    ENDIF 
    409      
    410     IF ( do_irrigation ) THEN  
    411        CALL xios_orchidee_send_field("irrigmap",irrigated) 
    412         
    413        IF (printlev >= 3) WRITE(numout,*) 'Verification : range of irrigated : ', MINVAL(irrigated), MAXVAL(irrigated)  
    414        IF ( .NOT. almaoutput ) THEN 
    415           CALL histwrite_p(hist_id, 'irrigmap', 1, irrigated, nbpt, index) 
    416        ELSE 
    417           CALL histwrite_p(hist_id, 'IrrigationMap', 1, irrigated, nbpt, index) 
    418        ENDIF 
    419        IF ( hist2_id > 0 ) THEN 
    420           IF ( .NOT. almaoutput ) THEN 
    421              CALL histwrite_p(hist2_id, 'irrigmap', 1, irrigated, nbpt, index) 
    422           ELSE 
    423              CALL histwrite_p(hist2_id, 'IrrigationMap', 1, irrigated, nbpt, index) 
    424           ENDIF 
    425        ENDIF 
     428 
     429    IF ( do_irrigation ) THEN 
     430 
     431       IF (printlev >= 3) WRITE(numout,*) 'Verification : range of irrigated : ', MINVAL(irrigated), MAXVAL(irrigated) 
     432       IF (printlev >= 3) WRITE(numout,*) 'Verification : range of irrig_frac : ', MINVAL(irrig_frac), MAXVAL(irrig_frac) 
    426433    ENDIF 
    427      
    428     IF ( do_floodplains ) THEN 
    429        CALL xios_orchidee_send_field("floodmap",floodplains) 
    430         
    431        IF (printlev>=3) WRITE(numout,*) 'Verification : range of floodplains : ', MINVAL(floodplains), MAXVAL(floodplains)  
    432        IF ( .NOT. almaoutput ) THEN 
    433           CALL histwrite_p(hist_id, 'floodmap', 1, floodplains, nbpt, index) 
    434        ELSE 
    435           CALL histwrite_p(hist_id, 'FloodplainsMap', 1, floodplains, nbpt, index) 
    436        ENDIF 
    437        IF ( hist2_id > 0 ) THEN 
    438           IF ( .NOT. almaoutput ) THEN 
    439              CALL histwrite_p(hist2_id, 'floodmap', 1, floodplains, nbpt, index) 
    440           ELSE 
    441              CALL histwrite_p(hist2_id, 'FloodplainsMap', 1, floodplains, nbpt, index) 
    442           ENDIF 
    443        ENDIF 
    444     ENDIF 
    445      
     434 
    446435    IF ( doswamps ) THEN 
    447436       CALL xios_orchidee_send_field("swampmap",swamp) 
    … …  
    558547       & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, & 
    559548       & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, & 
    560        & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id) 
     549       & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id,& 
     550       soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw) ! 
    561551 
    562552    IMPLICIT NONE 
    … …  
    584574    REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile 
    585575    REAL(r_std), INTENT(in)        :: reinf_slope(nbpt)    !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1) 
     576    REAL(r_std), INTENT(in)        :: root_deficit(nbpt)   !! soil water deficit 
     577    REAL(r_std), INTENT(in)        :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless) 
     578    REAL(r_std), INTENT(in)        :: irrig_frac(nbpt)     !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE. 
     579    REAL(r_std), INTENT(in)        :: irrigated_next (nbpt)!! Dynamic irrig. area, calculated in slowproc and passed to routing 
     580    REAL(r_std), INTENT(in)        :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac 
    586581 
    587582    !! 0.2 Output variables 
    … …  
    670665    ! the growing seasons we will give more weight to these areas. 
    671666    ! 
    672     DO jv=2,nvm 
    673        DO ig=1,nbpt 
    674           humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing 
    675           vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing 
    676        ENDDO 
    677     ENDDO 
     667    ! New irrigation scheme uses mean of vegtot with jv 1 to nvm 
     668    ! Old scheme keeps jv 2 to nvm, even if possibly an error 
     669    IF ( .NOT. old_irrig_scheme ) THEN 
     670      DO jv=1,nvm 
     671         DO ig=1,nbpt 
     672            humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing 
     673            vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing 
     674         ENDDO 
     675      ENDDO 
     676    ELSE 
     677      DO jv=2,nvm 
     678         DO ig=1,nbpt 
     679            humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing 
     680            vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing 
     681         ENDDO 
     682      ENDDO 
     683    ENDIF 
     684    ! Here updates irrigmap to irrigated_next from slowproc, every timestep 
     685    !irrigated_next is updated in slowproc when time comes 
     686    !irrig_frac was also updated in slowproc, here used as input variable 
     687 
     688    IF ( do_irrigation .AND. irrig_map_dynamic_flag  ) THEN 
     689        irrigated(:) = irrigated_next(:) 
     690    ENDIF 
    678691    ! 
    679692    time_counter = time_counter + dt_sechiba  
    … …  
    685698       !! Computes the transport of water in the various reservoirs 
    686699       ! 
    687        CALL routing_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, & 
     700       CALL routing_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, root_deficit, soiltile, & 
    688701            & vegtot_mean, totnobio_mean, transpot_mean, precip_mean, humrel_mean, k_litt_mean, floodtemp, reinf_slope, & 
    689             & lakeinflow_mean, returnflow_mean, reinfiltration_mean, irrigation_mean, riverflow_mean, & 
     702            & lakeinflow_mean, returnflow_mean, reinfiltration_mean, irrig_frac, irrigation_mean, irrigdeficit_mean, & 
     703            & irrigadduct_mean, irrig_gw_source_mean, & ! 
     704            & irrig_fast_source_mean, irrig_str_source_mean, riverflow_mean, & ! 
    690705            & coastalflow_mean, hydrographs, slowflow_diag, flood_frac, flood_res, & 
    691             & netflow_stream_diag, netflow_fast_diag, netflow_slow_diag) 
     706            & netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, fraction_aeirrig_sw) 
    692707       ! 
    693708       !! Responsible for storing the water in lakes 
    … …  
    718733       reinfiltration_mean(:) = reinfiltration_mean(:)/dt_routing*dt_sechiba 
    719734       irrigation_mean(:) = irrigation_mean(:)/dt_routing*dt_sechiba 
     735       irrigdeficit_mean(:) = irrigdeficit_mean/dt_routing*dt_sechiba ! 
     736       irrigadduct_mean(:) = irrigadduct_mean(:)/dt_routing*dt_sechiba! 
     737       irrig_gw_source_mean(:) = irrig_gw_source_mean(:)/dt_routing*dt_sechiba ! 
     738       irrig_fast_source_mean(:) = irrig_fast_source_mean(:)/dt_routing*dt_sechiba ! 
     739       irrig_str_source_mean(:) = irrig_str_source_mean(:)/dt_routing*dt_sechiba ! 
    720740       irrig_netereq(:) = irrig_netereq(:)/dt_routing*dt_sechiba 
    721741        
    … …  
    766786    ! Water fluxes converted from kg/m^2/dt_sechiba into kg/m^2/s  
    767787    CALL xios_orchidee_send_field("irrigation",irrigation/dt_sechiba) 
     788    CALL xios_orchidee_send_field("irrig_deficit",irrigdeficit_mean/dt_sechiba)! 
     789    CALL xios_orchidee_send_field("irrig_adduct",irrigadduct_mean/dt_sechiba)! 
     790    CALL xios_orchidee_send_field("irrig_gw_source",irrig_gw_source_mean/dt_sechiba) ! 
     791    CALL xios_orchidee_send_field("irrig_fast_source",irrig_fast_source_mean/dt_sechiba) ! 
     792    CALL xios_orchidee_send_field("irrig_str_source",irrig_str_source_mean/dt_sechiba) ! 
    768793    CALL xios_orchidee_send_field("netirrig",irrig_netereq/dt_sechiba) 
    769794    CALL xios_orchidee_send_field("riversret",returnflow/dt_sechiba) 
    … …  
    879904     
    880905    !! 0.2 Local variables 
    881     REAL(r_std), DIMENSION(1)      :: tmp_day               
     906    REAL(r_std), DIMENSION(1)      :: tmp_day 
    882907 
    883908!_ ================================================================================================================================ 
    … …  
    931956 
    932957    IF ( do_irrigation ) THEN 
    933        CALL restput_p (rest_id, 'irrigated', nbp_glo, 1, 1, kjit, irrigated, 'scatter',  nbp_glo, index_g) 
     958       !CALL restput_p (rest_id, 'irrigated', nbp_glo, 1, 1, kjit, irrigated, 'scatter',  nbp_glo, index_g) 
    934959       CALL restput_p (rest_id, 'irrigation', nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter',  nbp_glo, index_g) 
     960       CALL restput_p (rest_id, 'irrigdeficit', nbp_glo, 1, 1, kjit, irrigdeficit_mean, 'scatter',  nbp_glo, index_g) 
     961       CALL restput_p (rest_id, 'irrigadduct', nbp_glo, 1, 1, kjit, irrigadduct_mean, 'scatter',  nbp_glo, index_g) 
     962       CALL restput_p (rest_id, 'irrig_gw_source', nbp_glo, 1, 1, kjit, irrig_gw_source_mean, 'scatter',  nbp_glo, index_g) ! 
     963       CALL restput_p (rest_id, 'irrig_fast_source', nbp_glo, 1, 1, kjit, irrig_fast_source_mean, 'scatter',  nbp_glo, index_g) ! 
     964       CALL restput_p (rest_id, 'irrig_str_source', nbp_glo, 1, 1, kjit, irrig_str_source_mean, 'scatter',  nbp_glo, index_g) ! 
    935965    ENDIF 
    936966 
    … …  
    10221052    ! 
    10231053    !Config Key   = DO_FLOODINFILT 
    1024     !Config Desc  = Should floodplains reinfiltrate into the soil  
     1054    !Config Desc  = Should floodplains reinfiltrate into the soil 
    10251055    !Config If    = RIVER_ROUTING 
    10261056    !Config Def   = n 
    10271057    !Config Help  = This parameters allows the user to ask the model 
    1028     !Config         to take into account the flood plains reinfiltration  
    1029     !Config         into the soil moisture. It then can go  
     1058    !Config         to take into account the flood plains reinfiltration 
     1059    !Config         into the soil moisture. It then can go 
    10301060    !Config         back to the slow and fast reservoirs 
    10311061    !Config Units = [FLAG] 
    … …  
    10391069    !Config Def   = n 
    10401070    !Config Help  = This parameters allows the user to ask the model 
    1041     !Config         to take into account the swamps and return  
    1042     !Config         the water into the bottom of the soil. It then can go  
    1043     !Config         back to the atmopshere. This tried to simulate  
     1071    !Config         to take into account the swamps and return 
     1072    !Config         the water into the bottom of the soil. It then can go 
     1073    !Config         back to the atmopshere. This tried to simulate 
    10441074    !Config         internal deltas of rivers. 
    10451075    !Config Units = [FLAG] 
    … …  
    10491079    ! 
    10501080    !Config Key   = DO_PONDS 
    1051     !Config Desc  = Should we include ponds  
     1081    !Config Desc  = Should we include ponds 
    10521082    !Config If    = RIVER_ROUTING 
    10531083    !Config Def   = n 
    10541084    !Config Help  = This parameters allows the user to ask the model 
    1055     !Config         to take into account the ponds and return  
    1056     !Config         the water into the soil moisture. It then can go  
    1057     !Config         back to the atmopshere. This tried to simulate  
     1085    !Config         to take into account the ponds and return 
     1086    !Config         the water into the soil moisture. It then can go 
     1087    !Config         back to the atmopshere. This tried to simulate 
    10581088    !Config         little ponds especially in West Africa. 
    10591089    !Config Units = [FLAG] 
    … …  
    10641094 
    10651095    !Config Key   = SLOW_TCST 
    1066     !Config Desc  = Time constant for the slow reservoir  
    1067     !Config If    = RIVER_ROUTING  
    1068     !Config Def   = 25.0  
    1069     !Config Help  = This parameters allows the user to fix the  
     1096    !Config Desc  = Time constant for the slow reservoir 
     1097    !Config If    = RIVER_ROUTING 
     1098    !Config Def   = 25.0 
     1099    !Config Help  = This parameters allows the user to fix the 
    10701100    !Config         time constant (in days) of the slow reservoir 
    1071     !Config         in order to get better river flows for  
     1101    !Config         in order to get better river flows for 
    10721102    !Config         particular regions. 
    10731103    !Config Units = [days] 
    … …  
    10781108!> during the 1 degree NCEP Corrected by Cru (NCC) resolution simulations (Ngo-Duc et al., 2005, Ngo-Duc et al., 2006) and 
    10791109!> generalized for all the basins of the world. The "slow reservoir" and the "fast reservoir" 
    1080 !> have the highest value in order to simulate the groundwater.  
     1110!> have the highest value in order to simulate the groundwater. 
    10811111!> The "stream reservoir", which represents all the water of the stream, has the lowest value. 
    10821112!> Those figures are the same for all the basins of the world. 
    … …  
    10871117    ! 
    10881118    !Config Key   = FAST_TCST 
    1089     !Config Desc  = Time constant for the fast reservoir  
    1090     !Config If    = RIVER_ROUTING  
    1091     !Config Def   = 3.0  
    1092     !Config Help  = This parameters allows the user to fix the  
     1119    !Config Desc  = Time constant for the fast reservoir 
     1120    !Config If    = RIVER_ROUTING 
     1121    !Config Def   = 3.0 
     1122    !Config Help  = This parameters allows the user to fix the 
    10931123    !Config         time constant (in days) of the fast reservoir 
    1094     !Config         in order to get better river flows for  
     1124    !Config         in order to get better river flows for 
    10951125    !Config         particular regions. 
    10961126    !Config Units = [days] 
    10971127    CALL getin_p('FAST_TCST', fast_tcst) 
    1098      
     1128 
    10991129    !Config Key   = STREAM_TCST 
    1100     !Config Desc  = Time constant for the stream reservoir  
     1130    !Config Desc  = Time constant for the stream reservoir 
    11011131    !Config If    = RIVER_ROUTING 
    11021132    !Config Def   = 0.24 
    1103     !Config Help  = This parameters allows the user to fix the  
     1133    !Config Help  = This parameters allows the user to fix the 
    11041134    !Config         time constant (in days) of the stream reservoir 
    1105     !Config         in order to get better river flows for  
     1135    !Config         in order to get better river flows for 
    11061136    !Config         particular regions. 
    11071137    !Config Units = [days] 
    11081138    CALL getin_p('STREAM_TCST', stream_tcst) 
    1109      
     1139 
    11101140    !Config Key   = FLOOD_TCST 
    1111     !Config Desc  = Time constant for the flood reservoir  
     1141    !Config Desc  = Time constant for the flood reservoir 
    11121142    !Config If    = RIVER_ROUTING 
    11131143    !Config Def   = 4.0 
    1114     !Config Help  = This parameters allows the user to fix the  
     1144    !Config Help  = This parameters allows the user to fix the 
    11151145    !Config         time constant (in days) of the flood reservoir 
    1116     !Config         in order to get better river flows for  
     1146    !Config         in order to get better river flows for 
    11171147    !Config         particular regions. 
    11181148    !Config Units = [days] 
    11191149    CALL getin_p('FLOOD_TCST', flood_tcst) 
    1120      
     1150 
    11211151    !Config Key   = SWAMP_CST 
    1122     !Config Desc  = Fraction of the river that flows back to swamps  
     1152    !Config Desc  = Fraction of the river that flows back to swamps 
    11231153    !Config If    = RIVER_ROUTING 
    11241154    !Config Def   = 0.2 
    1125     !Config Help  = This parameters allows the user to fix the  
     1155    !Config Help  = This parameters allows the user to fix the 
    11261156    !Config         fraction of the river transport 
    11271157    !Config         that flows to swamps 
    11281158    !Config Units = [-] 
    11291159    CALL getin_p('SWAMP_CST', swamp_cst) 
    1130      
     1160 
    11311161    !Config Key   = FLOOD_BETA 
    1132     !Config Desc  = Parameter to fix the shape of the floodplain   
     1162    !Config Desc  = Parameter to fix the shape of the floodplain 
    11331163    !Config If    = RIVER_ROUTING 
    11341164    !Config Def   = 2.0 
    11351165    !Config Help  = Parameter to fix the shape of the floodplain 
    11361166    !Config         (>1 for convex edges, <1 for concave edges) 
    1137     !Config Units = [-]  
     1167    !Config Units = [-] 
    11381168    CALL getin_p("FLOOD_BETA", beta) 
    11391169    ! 
    … …  
    11421172    !Config If    = RIVER_ROUTING 
    11431173    !Config Def   = 0.5 
    1144     !Config Help  =  
    1145     !Config Units = [-]  
    1146     CALL getin_p("POND_BETAP", betap)     
     1174    !Config Help  = 
     1175    !Config Units = [-] 
     1176    CALL getin_p("POND_BETAP", betap) 
    11471177    ! 
    11481178    !Config Key   = FLOOD_CRI 
    … …  
    11501180    !Config If    = DO_FLOODPLAINS or DO_PONDS 
    11511181    !Config Def   = 2000. 
    1152     !Config Help  =  
    1153     !Config Units = [mm]  
     1182    !Config Help  = 
     1183    !Config Units = [mm] 
    11541184    CALL getin_p("FLOOD_CRI", floodcri) 
    11551185    ! 
    … …  
    11581188    !Config If    = DO_FLOODPLAINS or DO_PONDS 
    11591189    !Config Def   = 2000. 
    1160     !Config Help  =  
    1161     !Config Units = [mm]  
     1190    !Config Help  = 
     1191    !Config Units = [mm] 
    11621192    CALL getin_p("POND_CRI", pondcri) 
    11631193 
    … …  
    11661196    !Config If    = RIVER_ROUTING 
    11671197    !Config Def   = 7000 
    1168     !Config Help  =  
    1169     !Config Units = [kg/m2(routing area)]  
     1198    !Config Help  = 
     1199    !Config Units = [kg/m2(routing area)] 
    11701200    max_lake_reservoir = 7000 
    11711201    CALL getin_p("MAX_LAKE_RESERVOIR", max_lake_reservoir) 
    … …  
    12051235       ELSE 
    12061236          ! Take the value from restart file 
    1207           time_counter = tmp_day(1)  
     1237          time_counter = tmp_day(1) 
    12081238       ENDIF 
    12091239    ENDIF 
    12101240    CALL bcast(time_counter) 
    12111241 
    1212      
     1242 
    12131243    ALLOCATE (routing_area_loc(nbpt,nbasmax), stat=ier) 
    12141244    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for routing_area_loc','','') 
    … …  
    13651395    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_height, "gather", nbp_glo, index_g) 
    13661396    CALL setvar_p (flood_height, val_exp, 'NO_KEYWORD', zero) 
    1367      
     1397 
    13681398    ALLOCATE (pond_frac(nbpt), stat=ier) 
    13691399    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for pond_frac','','') 
    … …  
    13731403    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_frac, "gather", nbp_glo, index_g) 
    13741404    CALL setvar_p (pond_frac, val_exp, 'NO_KEYWORD', zero) 
    1375      
     1405 
    13761406    var_name = 'flood_frac' 
    13771407    CALL ioconf_setatt_p('UNITS', '-') 
    … …  
    13791409    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_frac, "gather", nbp_glo, index_g) 
    13801410    CALL setvar_p (flood_frac, val_exp, 'NO_KEYWORD', zero) 
    1381      
     1411 
    13821412    var_name = 'flood_res' 
    13831413    CALL ioconf_setatt_p('UNITS','mm') 
    … …  
    13861416    CALL setvar_p (flood_res, val_exp, 'NO_KEYWORD', zero) 
    13871417!    flood_res = zero 
    1388      
     1418 
    13891419    ALLOCATE (lake_reservoir(nbpt), stat=ier) 
    13901420    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for lake_reservoir','','') 
    … …  
    13941424    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lake_reservoir, "gather", nbp_glo, index_g) 
    13951425    CALL setvar_p (lake_reservoir, val_exp, 'NO_KEYWORD', zero) 
    1396      
     1426 
    13971427    ALLOCATE (pond_reservoir(nbpt), stat=ier) 
    13981428    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for pond_reservoir','','') 
    … …  
    14041434    ! 
    14051435    ! Map of irrigated areas 
    1406     ! 
    1407     IF ( do_irrigation ) THEN 
    1408        ALLOCATE (irrigated(nbpt), stat=ier) 
    1409        IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigated','','') 
    1410        var_name = 'irrigated' 
    1411        CALL ioconf_setatt_p('UNITS', 'm^2') 
    1412        CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area') 
    1413        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigated, "gather", nbp_glo, index_g) 
    1414        CALL setvar_p (irrigated, val_exp, 'NO_KEYWORD', undef_sechiba) 
    1415     ENDIF 
    1416      
     1436    ! irrigated equal to irrigated_next from slowproc. Here, we initialize to zero. 
     1437    ! Values from slowproc given in routing_main 
     1438    ALLOCATE (irrigated(nbpt), stat=ier) 
     1439    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigated','','') 
     1440    var_name = 'irrigated' 
     1441    CALL ioconf_setatt_p('UNITS', 'm^2') 
     1442    CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area') 
     1443    irrigated(:) = zero 
     1444 
    14171445    IF ( do_floodplains ) THEN 
    14181446       ALLOCATE (floodplains(nbpt), stat=ier) 
    … …  
    14431471    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g) 
    14441472    CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero) 
    1445      
     1473 
    14461474    ALLOCATE (returnflow_mean(nbpt), stat=ier) 
    14471475    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for returnflow_mean','','') 
    … …  
    14521480    CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', zero) 
    14531481    returnflow(:) = returnflow_mean(:) 
    1454      
     1482 
    14551483    ALLOCATE (reinfiltration_mean(nbpt), stat=ier) 
    14561484    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for reinfiltration_mean','','') 
    … …  
    14611489    CALL setvar_p (reinfiltration_mean, val_exp, 'NO_KEYWORD', zero) 
    14621490    reinfiltration(:) = reinfiltration_mean(:) 
    1463      
     1491 
    14641492    ALLOCATE (irrigation_mean(nbpt), stat=ier) 
    14651493    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigation_mean','','') 
     1494    ALLOCATE (irrigdeficit_mean(nbpt), stat=ier) 
     1495    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigdeficit_mean','','') 
     1496    ALLOCATE (irrigadduct_mean(nbpt), stat=ier) 
     1497    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigadduct_mean','','') 
    14661498    ALLOCATE (irrig_netereq(nbpt), stat=ier) 
    14671499    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_netereq','','') 
     1500    ALLOCATE (irrig_gw_source_mean(nbpt), stat=ier) ! 
     1501    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_gw_source_mean','','')! 
     1502    ALLOCATE (irrig_fast_source_mean(nbpt), stat=ier) ! 
     1503    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_fast_source_mean','','')! 
     1504    ALLOCATE (irrig_str_source_mean(nbpt), stat=ier) ! 
     1505    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_str_source_mean','','')! 
    14681506    irrig_netereq(:) = zero 
    1469      
     1507 
    14701508    IF ( do_irrigation ) THEN 
    14711509       var_name = 'irrigation' 
    … …  
    14741512       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g) 
    14751513       CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero) 
     1514       var_name = 'irrigdeficit' 
     1515       CALL ioconf_setatt_p('UNITS', 'Kg/dt') 
     1516       CALL ioconf_setatt_p('LONG_NAME','Irrigation deficit') 
     1517       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigdeficit_mean, "gather", nbp_glo, index_g) 
     1518       CALL setvar_p (irrigdeficit_mean, val_exp, 'NO_KEYWORD', zero) 
     1519       var_name = 'irrigadduct' 
     1520       CALL ioconf_setatt_p('UNITS', 'Kg/dt') 
     1521       CALL ioconf_setatt_p('LONG_NAME','Artificial irrigation flux from adduction') 
     1522       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigadduct_mean, "gather", nbp_glo, index_g) 
     1523       CALL setvar_p (irrigadduct_mean, val_exp, 'NO_KEYWORD', zero) 
     1524       var_name = 'irrig_gw_source'! 
     1525       CALL ioconf_setatt_p('UNITS', 'Kg/dt') 
     1526       CALL ioconf_setatt_p('LONG_NAME','Irrigation from GW reservoir') 
     1527       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrig_gw_source_mean, "gather", nbp_glo, index_g) 
     1528       CALL setvar_p (irrig_gw_source_mean, val_exp, 'NO_KEYWORD', zero) 
     1529       var_name = 'irrig_fast_source'! 
     1530       CALL ioconf_setatt_p('UNITS', 'Kg/dt') 
     1531       CALL ioconf_setatt_p('LONG_NAME','Irrigation from fast reservoir') 
     1532       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrig_fast_source_mean, "gather", nbp_glo, index_g) 
     1533       CALL setvar_p (irrig_fast_source_mean, val_exp, 'NO_KEYWORD', zero) 
     1534       var_name = 'irrig_str_source'! 
     1535       CALL ioconf_setatt_p('UNITS', 'Kg/dt') 
     1536       CALL ioconf_setatt_p('LONG_NAME','Irrigation from stream reservoir') 
     1537       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrig_str_source_mean, "gather", nbp_glo, index_g) 
     1538       CALL setvar_p (irrig_str_source_mean, val_exp, 'NO_KEYWORD', zero) 
    14761539    ELSE 
    14771540       irrigation_mean(:) = zero 
     1541       irrig_gw_source_mean(:) = zero 
     1542       irrig_fast_source_mean(:) = zero 
     1543       irrig_str_source_mean(:) = zero 
     1544       irrigdeficit_mean(:) = zero 
     1545       irrigadduct_mean(:) = zero 
    14781546    ENDIF 
    1479     irrigation(:) = irrigation_mean(:)  
    1480      
     1547    irrigation(:) = irrigation_mean(:) 
     1548 
    14811549    ALLOCATE (riverflow_mean(nbpt), stat=ier) 
    14821550    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for riverflow_mean','','') 
    … …  
    14871555    CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero) 
    14881556    riverflow(:) = riverflow_mean(:) 
    1489      
     1557 
    14901558    ALLOCATE (coastalflow_mean(nbpt), stat=ier) 
    14911559    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for coastalflow_mean','','') 
    … …  
    14961564    CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero) 
    14971565    coastalflow(:) = coastalflow_mean(:) 
    1498      
     1566 
    14991567    ! Locate it at the 2m level 
    15001568    ipn = MINLOC(ABS(diaglev-2)) 
    … …  
    15031571    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for floodtemp','','') 
    15041572    floodtemp(:) = stempdiag(:,floodtemp_lev) 
    1505      
     1573 
    15061574    ALLOCATE(hydrographs(nbpt), stat=ier) 
    15071575    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for hydrographs','','') 
    … …  
    15111579    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs, "gather", nbp_glo, index_g) 
    15121580    CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', zero) 
    1513   
     1581 
    15141582    ALLOCATE(slowflow_diag(nbpt), stat=ier) 
    15151583    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for slowflow_diag','','') 
    … …  
    15261594         & pond_diag(nbpt), lake_diag(nbpt), stat=ier) 
    15271595    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for fast_diag,..','','') 
    1528      
     1596 
    15291597    fast_diag(:) = zero 
    15301598    slow_diag(:) = zero 
    … …  
    15331601    pond_diag(:) = zero 
    15341602    lake_diag(:) = zero 
    1535      
     1603 
    15361604    DO ig=1,nbpt 
    15371605       totarea = zero 
    … …  
    15661634    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodout_mean, "gather", nbp_glo, index_g) 
    15671635    CALL setvar_p (floodout_mean, val_exp, 'NO_KEYWORD', zero) 
    1568      
     1636 
    15691637    ALLOCATE (runoff_mean(nbpt), stat=ier) 
    15701638    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for runoff_mean','','') 
    … …  
    15741642    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g) 
    15751643    CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero) 
    1576      
     1644 
    15771645    ALLOCATE(drainage_mean(nbpt), stat=ier) 
    15781646    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for drainage_mean','','') 
    … …  
    15821650    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g) 
    15831651    CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero) 
    1584      
     1652 
    15851653    ALLOCATE(transpot_mean(nbpt), stat=ier) 
    15861654    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for transpot_mean','','') 
    … …  
    15981666    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., precip_mean, "gather", nbp_glo, index_g) 
    15991667    CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', zero) 
    1600      
     1668 
    16011669    ALLOCATE(humrel_mean(nbpt), stat=ier) 
    16021670    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for humrel_mean','','') 
    … …  
    16061674    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., humrel_mean, "gather", nbp_glo, index_g) 
    16071675    CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', un) 
    1608      
     1676 
    16091677    ALLOCATE(k_litt_mean(nbpt), stat=ier) 
    16101678    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for k_litt_mean','','') 
    … …  
    16141682    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., k_litt_mean, "gather", nbp_glo, index_g) 
    16151683    CALL setvar_p (k_litt_mean, val_exp, 'NO_KEYWORD', zero) 
    1616      
     1684 
    16171685    ALLOCATE(totnobio_mean(nbpt), stat=ier) 
    16181686    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for totnobio_mean','','') 
    … …  
    16221690    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., totnobio_mean, "gather", nbp_glo, index_g) 
    16231691    CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', zero) 
    1624      
     1692 
    16251693    ALLOCATE(vegtot_mean(nbpt), stat=ier) 
    16261694    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for vegtot_mean','','') 
    … …  
    16931761    IF (ALLOCATED(hydrodiag_loc)) DEALLOCATE(hydrodiag_loc) 
    16941762    IF (ALLOCATED(hydrodiag_glo)) DEALLOCATE(hydrodiag_glo) 
    1695     IF (ALLOCATED(hydroupbasin_loc)) DEALLOCATE(hydroupbasin_loc)     
     1763    IF (ALLOCATED(hydroupbasin_loc)) DEALLOCATE(hydroupbasin_loc) 
    16961764    IF (ALLOCATED(hydroupbasin_glo)) DEALLOCATE(hydroupbasin_glo) 
    16971765    IF (ALLOCATED(hydrographs)) DEALLOCATE(hydrographs) 
    16981766    IF (ALLOCATED(slowflow_diag)) DEALLOCATE(slowflow_diag) 
    16991767    IF (ALLOCATED(irrigation_mean)) DEALLOCATE(irrigation_mean) 
     1768    IF (ALLOCATED(irrigdeficit_mean)) DEALLOCATE(irrigdeficit_mean)! 
     1769    IF (ALLOCATED(irrigadduct_mean)) DEALLOCATE(irrigadduct_mean)! 
     1770    IF (ALLOCATED(irrig_gw_source_mean)) DEALLOCATE(irrig_gw_source_mean) ! 
     1771    IF (ALLOCATED(irrig_fast_source_mean)) DEALLOCATE(irrig_fast_source_mean) ! 
     1772    IF (ALLOCATED(irrig_str_source_mean)) DEALLOCATE(irrig_str_source_mean) ! 
    17001773    IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated) 
    17011774    IF (ALLOCATED(floodplains)) DEALLOCATE(floodplains) 
    … …  
    17181791!! 
    17191792!! DESCRIPTION (definitions, functional, design, flags) : 
    1720 !! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an  
    1721 !! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes  
    1722 !! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of  
     1793!! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an 
     1794!! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes 
     1795!! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of 
    17231796!! the up-stream basin and adding it to the down-stream one. 
    1724 !! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once  
     1797!! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once 
    17251798!! all reservoirs are updated we deal with irrigation. The final step is to compute diagnostic fluxes. Among them 
    17261799!! the hydrographs of the largest rivers we have chosen to monitor. 
    … …  
    17781851!_ ================================================================================================================================ 
    17791852 
    1780   SUBROUTINE routing_flow(nbpt, dt_routing, lalo, floodout, runoff, drainage, & 
     1853  SUBROUTINE routing_flow(nbpt, dt_routing, lalo, floodout, runoff, drainage, root_deficit, soiltile, & 
    17811854       &                  vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, reinf_slope, & 
    1782        &                  lakeinflow, returnflow, reinfiltration, irrigation, riverflow, & 
     1855       &                  lakeinflow, returnflow, reinfiltration, irrig_frac, irrigation, irrigdeficit, & 
     1856       &                  irrigadduct, irrig_gw_source, & 
     1857       &                  irrig_fast_source, irrig_str_source, riverflow, & ! 
    17831858       &                  coastalflow, hydrographs, slowflow_diag, flood_frac, flood_res, & 
    1784                           netflow_stream_diag, netflow_fast_diag, netflow_slow_diag) 
     1859                          netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, fraction_aeirrig_sw) ! 
    17851860    ! 
    17861861    IMPLICIT NONE 
    … …  
    17931868    REAL(r_std), INTENT(in)                      :: floodout(nbpt)            !! Grid-point flow out of floodplains (kg/m^2/dt) 
    17941869    REAL(r_std), INTENT(in)                      :: drainage(nbpt)            !! Grid-point drainage (kg/m^2/dt) 
     1870    REAL(r_std), INTENT(in)                      :: root_deficit(nbpt)        !! soil water deficit 
    17951871    REAL(r_std), INTENT(in)                      :: vegtot(nbpt)              !! Potentially vegetated fraction (unitless;0-1) 
    17961872    REAL(r_std), INTENT(in)                      :: totnobio(nbpt)            !! Other areas which can not have vegetation 
    … …  
    18021878    REAL(r_std), INTENT(in)                      :: reinf_slope(nbpt)         !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1) 
    18031879    REAL(r_std), INTENT(out)                     :: lakeinflow(nbpt)          !! Water inflow to the lakes (kg/dt) 
     1880    REAL(r_std), INTENT(in)                      :: soiltile(nbpt,nstm)       !! Fraction of each soil tile within vegtot (0-1, unitless) 
     1881    REAL(r_std), INTENT(in)                      :: irrig_frac(nbpt)          !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE. 
     1882    REAL(r_std), INTENT(in)                      :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac 
     1883 
    18041884    ! 
    18051885!! OUTPUT VARIABLES 
    … …  
    18081888    REAL(r_std), INTENT(out)                     :: reinfiltration(nbpt)      !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) 
    18091889    REAL(r_std), INTENT(out)                     :: irrigation(nbpt)          !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt_routing) 
     1890    REAL(r_std), INTENT(out)                     :: irrigdeficit(nbpt)          !! Irrigation deficit. Difference btw irrig. demand and irrigation 
     1891    REAL(r_std), INTENT(out)                     :: irrigadduct(nbpt)          !! Irrigation from adducted water. Included in irrigation 
     1892 
    18101893    REAL(r_std), INTENT(out)                     :: riverflow(nbpt)           !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt_routing) 
    18111894    REAL(r_std), INTENT(out)                     :: coastalflow(nbpt)         !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt_routing) 
    … …  
    18181901    REAL(r_std), INTENT(out)                     :: netflow_fast_diag(nbpt)   !! Input - Output flow to fast reservoir 
    18191902    REAL(r_std), INTENT(out)                     :: netflow_slow_diag(nbpt)   !! Input - Output flow to slow reservoir 
     1903 
     1904    REAL(r_std), INTENT(out)                     :: irrig_gw_source(nbpt)     !! Irrigation from  groundwater per cell 
     1905    REAL(r_std), INTENT(out)                     :: irrig_fast_source(nbpt)     !! Irrigation from  fast reservoir per cell 
     1906    REAL(r_std), INTENT(out)                     :: irrig_str_source(nbpt)     !! Irrigation from  stramflow per cell 
    18201907    ! 
    18211908!! LOCAL VARIABLES 
    … …  
    18431930    REAL(r_std), DIMENSION(nbpt)                 :: totarea                   !! Total area of basin (m^2) 
    18441931    REAL(r_std), DIMENSION(nbpt)                 :: totflood                  !! Total amount of water in the floodplains reservoir (kg) 
    1845     REAL(r_std), DIMENSION(nbasmax)              :: pond_excessflow           !!  
     1932    REAL(r_std), DIMENSION(nbasmax)              :: pond_excessflow           !! 
    18461933    REAL(r_std)                                  :: flow                      !! Outflow computation for the reservoirs (kg/dt) 
    18471934    REAL(r_std)                                  :: floodindex                !! Fraction of grid box area inundated (unitless;0-1) 
    1848     REAL(r_std)                                  :: pondex                    !!  
     1935    REAL(r_std)                                  :: pondex                    !! 
    18491936    REAL(r_std)                                  :: flood_frac_pot            !! Total fraction of the grid box which is flooded at optimum repartition (unitless;0-1) 
    18501937    REAL(r_std)                                  :: stream_tot                !! Total water amount in the stream reservoirs (kg) 
    … …  
    18551942    REAL(r_std)                                  :: total_lake_overflow       !! Sum of lake_overflow over full grid (kg) 
    18561943    REAL(r_std), DIMENSION(8,nbasmax)            :: streams_around            !! Stream reservoirs of the neighboring grid boxes (kg) 
    1857     INTEGER(i_std), DIMENSION(8)                 :: igrd                      !!  
    1858     INTEGER(i_std), DIMENSION(2)                 :: ff                        !!  
    1859     INTEGER(i_std), DIMENSION(1)                 :: fi                        !!  
     1944    INTEGER(i_std), DIMENSION(8)                 :: igrd                      !! 
     1945    INTEGER(i_std), DIMENSION(2)                 :: ff                        !! 
     1946    INTEGER(i_std), DIMENSION(1)                 :: fi                        !! 
    18601947    INTEGER(i_std)                               :: ig, ib, ib2, ig2          !! Indices (unitless) 
    18611948    INTEGER(i_std)                               :: rtg, rtb, in              !! Indices (unitless) 
    1862     INTEGER(i_std)                               :: ier                       !! Error handling  
     1949    INTEGER(i_std)                               :: ier                       !! Error handling 
     1950    REAL(r_std), DIMENSION(nbpt)              :: Count_failure_slow        !! Counter times slow reserv. does not fit irrigation demand 
     1951    REAL(r_std), DIMENSION(nbpt)              :: Count_failure_fast        !! Counter times fast reserv. does not fit irrigation demand 
     1952    REAL(r_std), DIMENSION(nbpt)              :: Count_failure_stre        !! Counter times stream reserv. does not fit irrigation demand 
     1953    LOGICAL              :: IsFail_slow        !! Logical to ask if slow reserv. does not fit irrigation demand 
     1954    LOGICAL              :: IsFail_fast        !! Logical to ask if fast reserv. does not fit irrigation demand 
     1955    LOGICAL              :: IsFail_stre        !! Logical to ask if stream reserv. does not fit irrigation demand 
    18631956    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: fast_flow_g               !! Outflow from the fast reservoir (kg/dt) 
    18641957    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: slow_flow_g               !! Outflow from the slow reservoir (kg/dt) 
    … …  
    18721965    REAL(r_std), DIMENSION(nbpt, nbasmax)        :: netflow_slow              !! Input - Output flow to slow reservoir 
    18731966 
    1874  
     1967    REAL(r_std)                :: pcent_vol_irrig               !! Percentage of surface water for irrigation from total available water 
     1968    REAL(r_std)                :: slow_wdr_dummy,fast_wdr_dummy,stre_wdr_dummy !! Dummy variables, real abstraction in each reservoir for irrigation 
     1969    REAL(r_std)                :: pot_slow_wdr_dummy,pot_fast_wdr_dummy,pot_stre_wdr_dummy !! Dummy variables, potential abstraction in each reservoir for irrigatio 
    18751970    !! PARAMETERS 
    18761971    LOGICAL, PARAMETER                           :: check_reservoir = .FALSE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false) 
    … …  
    18861981    totarea(:) = zero 
    18871982    totflood(:) = zero 
     1983    irrig_gw_source(:) = zero ! 
     1984    irrig_fast_source(:) = zero ! 
     1985    irrig_str_source(:) = zero ! 
     1986    Count_failure_slow(:) = zero ! 
     1987    Count_failure_fast(:) = zero ! 
     1988    Count_failure_stre(:) = zero ! 
    18881989    ! 
    18891990    ! Compute all the fluxes 
    … …  
    18971998    ENDDO 
    18981999          ! 
    1899 !> The outflow fluxes from the three reservoirs are computed.  
     2000!> The outflow fluxes from the three reservoirs are computed. 
    19002001!> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume. 
    19012002!> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid 
    … …  
    19122013          IF ( route_tobasin(ig,ib) .GT. 0 ) THEN 
    19132014             ! 
    1914              ! Each of the fluxes is limited by the water in the reservoir and a small margin  
     2015             ! Each of the fluxes is limited by the water in the reservoir and a small margin 
    19152016             ! (min_reservoir) to avoid rounding errors. 
    19162017             ! 
    … …  
    19232024             slow_flow(ig,ib) = MAX(flow, zero) 
    19242025 
    1925              flow = MIN(stream_reservoir(ig,ib)/((topo_resid(ig,ib)/1000.)*stream_tcst* &  
     2026             flow = MIN(stream_reservoir(ig,ib)/((topo_resid(ig,ib)/1000.)*stream_tcst* & 
    19262027                  & MAX(un-SQRT(flood_frac_bas(ig,ib)),min_sechiba)*one_day/dt_routing),& 
    19272028                  & stream_reservoir(ig,ib)-min_sechiba) 
    … …  
    19442045             flow = MIN(floodout(ig)*totarea(ig)*pond_frac(ig)/flood_frac(ig), pond_reservoir(ig)+totflood(ig)) 
    19452046             pondex = MAX(flow - pond_reservoir(ig), zero) 
    1946              pond_reservoir(ig) = pond_reservoir(ig) - (flow - pondex)  
     2047             pond_reservoir(ig) = pond_reservoir(ig) - (flow - pondex) 
    19472048             ! 
    19482049             ! If demand was over reservoir size, we will take it out from floodplains 
    … …  
    19812082    !- 
    19822083    !- Computing the drainage and outflow from floodplains 
    1983 !> Drainage from floodplains is depending on a averaged conductivity (k_litt)  
     2084!> Drainage from floodplains is depending on a averaged conductivity (k_litt) 
    19842085!> for saturated infiltration in the 'litter' layer. Flood_drainage will be 
    19852086!> a component of the total reinfiltration that leaves the routing scheme. 
    … …  
    19972098          DO ib=1,nbasmax 
    19982099             DO ig=1,nbpt 
    1999                 flood_drainage(ig,ib) = zero  
     2100                flood_drainage(ig,ib) = zero 
    20002101             ENDDO 
    20012102          ENDDO 
    … …  
    20452146             pond_drainage(ig,ib) = MIN(pond_reservoir(ig)*routing_area(ig,ib)/totarea(ig), & 
    20462147                  & pond_frac(ig)*routing_area(ig,ib)*k_litt(ig)*dt_routing/one_day) 
    2047              fast_flow(ig,ib) = fast_flow(ig,ib) - pond_inflow(ig,ib)  
     2148             fast_flow(ig,ib) = fast_flow(ig,ib) - pond_inflow(ig,ib) 
    20482149          ENDDO 
    20492150       ENDDO 
    … …  
    20732174    ENDIF 
    20742175    IF (ier /= 0) CALL ipslerr_p(3,'routing_flow','Pb in allocate for fast_flow_g','','') 
    2075         
     2176 
    20762177    CALL gather(fast_flow,fast_flow_g) 
    20772178    CALL gather(slow_flow,slow_flow_g) 
    … …  
    20922193 
    20932194    DEALLOCATE( fast_flow_g, slow_flow_g, stream_flow_g ) 
    2094     
     2195 
    20952196    CALL scatter(transport_glo,transport) 
    20962197 
    … …  
    21122213       DO ib=1,nbasmax 
    21132214          DO ig=1,nbpt 
    2114              potflood(ig,ib) = transport(ig,ib)  
     2215             potflood(ig,ib) = transport(ig,ib) 
    21152216             ! 
    21162217             IF ( tobeflooded(ig) > 0. .AND. potflood(ig,ib) > 0. .AND. floodtemp(ig) > tp_00 ) THEN 
    … …  
    21232224                return_swamp(ig,ib) = swamp_cst * potflood(ig,ib) * floodindex 
    21242225                ! 
    2125                 tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib)  
     2226                tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib) 
    21262227                ! 
    21272228             ENDIF 
    … …  
    21362237          IF (floodplains(ig) .GT. min_sechiba .AND. floodtemp(ig) .GT. tp_00) THEN 
    21372238             DO ib=1,nbasmax 
    2138                 floods(ig,ib) = transport(ig,ib) - return_swamp(ig,ib)  
     2239                floods(ig,ib) = transport(ig,ib) - return_swamp(ig,ib) 
    21392240             ENDDO 
    21402241          ENDIF 
    … …  
    21622263          ! 
    21632264          flood_reservoir(ig,ib) = flood_reservoir(ig,ib) + floods(ig,ib) - & 
    2164                & flood_flow(ig,ib)  
     2265               & flood_flow(ig,ib) 
    21652266          ! 
    21662267          pond_reservoir(ig) = pond_reservoir(ig) + pond_inflow(ig,ib) - pond_drainage(ig,ib) 
    … …  
    21702271                WRITE(numout,*) "WARNING : negative flood reservoir at :", ig, ib, ". Problem is being corrected." 
    21712272                WRITE(numout,*) "flood_reservoir, floods, flood_flow : ", flood_reservoir(ig,ib), floods(ig,ib), & 
    2172                      & flood_flow(ig,ib)  
     2273                     & flood_flow(ig,ib) 
    21732274             ENDIF 
    21742275             stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + flood_reservoir(ig,ib) 
    … …  
    22292330             flood_frac_pot = (totflood(ig) / (totarea(ig)*floodcri/(beta+un)))**(beta/(beta+un)) 
    22302331             flood_frac(ig) = MIN(floodplains(ig) / totarea(ig), flood_frac_pot) 
    2231              ! Then we diagnose the fraction for each basin with the size of its flood_reservoir  
     2332             ! Then we diagnose the fraction for each basin with the size of its flood_reservoir 
    22322333             ! (flood_frac_bas may be > 1) 
    22332334             DO ib=1,nbasmax 
    … …  
    22382339             ENDDO 
    22392340             ! We diagnose the maximum height of floodplain 
    2240              flood_height(ig) = (beta/(beta+1))*floodcri*(flood_frac(ig))**(un/beta) + totflood(ig)/(totarea(ig)*flood_frac(ig))  
     2341             flood_height(ig) = (beta/(beta+1))*floodcri*(flood_frac(ig))**(un/beta) + totflood(ig)/(totarea(ig)*flood_frac(ig)) 
    22412342             ! And finally add the pond surface 
    2242              pond_frac(ig) = MIN(un-flood_frac(ig), ((betap+1)*pond_reservoir(ig) / (pondcri*totarea(ig)))**(betap/(betap+1)) )  
     2343             pond_frac(ig) = MIN(un-flood_frac(ig), ((betap+1)*pond_reservoir(ig) / (pondcri*totarea(ig)))**(betap/(betap+1)) ) 
    22432344             flood_frac(ig) = flood_frac(ig) + pond_frac(ig) 
    22442345             ! 
    … …  
    22652366          DO ig=1,nbpt 
    22662367             returnflow(ig) =  returnflow(ig) + return_swamp(ig,ib) 
    2267              reinfiltration(ig) =  reinfiltration(ig) + pond_drainage(ig,ib) + flood_drainage(ig,ib)  
     2368             reinfiltration(ig) =  reinfiltration(ig) + pond_drainage(ig,ib) + flood_drainage(ig,ib) 
    22682369          ENDDO 
    22692370       ENDDO 
    … …  
    23042405! 
    23052406    IF ( do_irrigation ) THEN 
    2306        DO ig=1,nbpt 
    2307           ! 
    2308           IF ((vegtot(ig) .GT. min_sechiba) .AND. (humrel(ig) .LT. un-min_sechiba) .AND. & 
    2309                & (runoff(ig) .LT. min_sechiba) ) THEN 
    2310               
    2311              irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, transpot_mean(ig) - & 
     2407      DO ig=1,nbpt 
     2408            !It enters to the new irrigation module only if there is an irrigated fraction, if not irrig_netereq = zero for that cell 
     2409            IF ((irrig_frac(ig) .GT. min_sechiba) .AND. .NOT. old_irrig_scheme ) THEN 
     2410 
     2411                  irrig_netereq(ig) = irrig_netereq(ig) + MIN( irrig_dosmax/3600*dt_routing, & 
     2412                                      root_deficit(ig) ) * soiltile(ig, irrig_st) * vegtot(ig) 
     2413                  ! By definition, irrig_dosmax is in kg/m2 of soil tile/hour,dividing by 3600(seconds/hour) * DT_ROUTING  ! 
     2414                  ! = kg/m2 of soil tile/(routing timestep) 
     2415                  ! irrig_netereq(kg/m2 of grid cell / routing timstep ) is equal to 
     2416                  ! root_deficit (kg/m2 of soil tile) * soiltile*vegtot (fraction of soil tile at cell level) = kg/m2 of grid cell 
     2417 
     2418                  IF (.NOT. irrigated_soiltile .AND. ( soiltile(ig,irrig_st) .GT. min_sechiba ) .AND. (vegtot(ig) .GT. min_sechiba) ) THEN 
     2419                      ! Irrigated_soiltile asks if there is an independent soil tile for irrigated crops. If not, 
     2420                      ! actual volume calculated for irrig_netereq assumed that the whole SOILTILE was irrigated, but in this case 
     2421                      ! just a fraction of the irrig_st (irrigated soil tile, by default = 3) is actually irrigated, 
     2422                      ! and irrig_netereq needs to be reduced by irrig_frac/( soiltile * vegtot ) (note that it is max = 1 thanks to irrig_frac calculation in l. 424) 
     2423                      ! Demand(ST3)*irrig_frac/(soiltile(3)*vegtot) = irrig_netereq_In_ST3, then 
     2424                      !irrig_netereq_In_ST3 * (soiltile(3)*vegtot) = irrig_netereq at grid scale = Demand(ST3)*irrig_frac. 
     2425                      irrig_netereq(ig) = irrig_netereq(ig) * irrig_frac(ig) / ( soiltile(ig,irrig_st) * vegtot(ig) ) 
     2426                      !irrig_netereq = kg/m2 of grid cell 
     2427 
     2428                  ENDIF 
     2429            !Old irrigation scheme as in tag 2.0 
     2430            ELSE IF((vegtot(ig) .GT. min_sechiba) .AND. (humrel(ig) .LT. un-min_sechiba) .AND. & 
     2431                    & (runoff(ig) .LT. min_sechiba) .AND.  old_irrig_scheme) THEN 
     2432 
     2433                  irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, transpot_mean(ig) - & 
    23122434                  & (precip(ig)+reinfiltration(ig)) ) 
    2313               
    2314           ENDIF 
    2315           ! 
    2316           DO ib=1,nbasmax 
    2317              IF ( routing_area(ig,ib) .GT. 0 ) THEN 
    2318               
    2319                 irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib) 
    2320  
    2321                 irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),& 
    2322                      &   stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) ) 
    2323                  
    2324                 slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + & 
    2325                      & MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib)))) 
    2326  
    2327                 fast_reservoir(ig,ib) = MAX( zero, & 
    2328                      &  fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib))) 
    2329  
    2330                 stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib) ) 
    2331  
    2332                 irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib) 
    2333  
    2334              ENDIF 
    2335           ENDDO 
     2435 
     2436            ENDIF 
     2437 
     2438            DO ib=1,nbasmax 
     2439                IF ( routing_area(ig,ib) .GT. 0 ) THEN 
     2440 
     2441                    IF (.NOT. old_irrig_scheme .AND. select_source_irrig) THEN 
     2442 
     2443                      ! For   irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir 
     2444                      ! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir 
     2445 
     2446                      ! The new priorization scheme takes into account irrig. infrastructur according to GMIA map 
     2447                      ! It also withdraw water according to availability, it means that it wont seek for all the water in then 
     2448                      ! stream reservoir, even if this one could respond to the demand by itself 
     2449 
     2450                      pot_slow_wdr_dummy = ( 1 - fraction_aeirrig_sw(ig)) * avail_reserve(3)*slow_reservoir(ig,ib) 
     2451                      pot_fast_wdr_dummy = fraction_aeirrig_sw(ig) * avail_reserve(2)*fast_reservoir(ig,ib) 
     2452                      pot_stre_wdr_dummy = fraction_aeirrig_sw(ig) * avail_reserve(1) * stream_reservoir(ig,ib) 
     2453                      pcent_vol_irrig = zero 
     2454                      IsFail_slow = .FALSE. ! 
     2455                      IsFail_fast = .FALSE. ! 
     2456                      IsFail_stre = .FALSE. ! 
     2457                      irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib) 
     2458 
     2459                      irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),& 
     2460                            pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy) 
     2461 
     2462                       !!   additional IF to calculate pcent_vol_irrig, in the case the total avail. 
     2463                       !! water is zero, I.E. when there is no water in surface and fraction_ae = 0, 
     2464                       !! so GW is not taken into account 
     2465                       !! Note on pcent_vol_irrig: It correspond to the fraction of available water in surface, 
     2466                       !! considering environmental needs and irrigation equipement by source from map. It controls 
     2467                       !! how the source of water withdrawl, especially when requirements < available water 
     2468 
     2469                      IF (  (pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy)  .GT. min_sechiba ) THEN 
     2470 
     2471                          pcent_vol_irrig = ( pot_stre_wdr_dummy + pot_fast_wdr_dummy ) / & 
     2472                                ( pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy ) 
     2473 
     2474                          !Irrig_actual set to zero, because there is no available water. 
     2475                          !Put to avoid negative values due to problems in the Min function 
     2476                          irrig_actual(ig,ib) = MAX(irrig_actual(ig,ib), zero) 
     2477                      !Already zero because pcent_vol_irrig initialized to zero 
     2478                      !Put here to readability but not necessary 
     2479                      !ELSE 
     2480                      !    pcent_vol_irrig = zero 
     2481 
     2482                      ENDIF 
     2483 
     2484                      !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the 
     2485                      !Volume used for irrigation that comes from GW 
     2486                      ! Idem for irrig_fast_source and irrig_str_source 
     2487 
     2488                      slow_wdr_dummy = slow_reservoir(ig,ib) 
     2489                      slow_reservoir(ig,ib) = MAX( (un - ( un - fraction_aeirrig_sw(ig) ) * avail_reserve(3) ) * & 
     2490                                              slow_reservoir(ig,ib), slow_reservoir(ig,ib) + & 
     2491                                              MIN( - irrig_actual(ig,ib) * (un - pcent_vol_irrig ), & 
     2492                                              avail_reserve(2) * fraction_aeirrig_sw(ig) * fast_reservoir(ig,ib) + & 
     2493                                              MIN(zero, avail_reserve(1) * fraction_aeirrig_sw(ig) * stream_reservoir(ig,ib)  - & 
     2494                                              pcent_vol_irrig * irrig_actual(ig,ib) ) ) ) 
     2495 
     2496                      slow_wdr_dummy = slow_wdr_dummy - slow_reservoir(ig,ib) 
     2497                      irrig_gw_source(ig) = irrig_gw_source(ig) + slow_wdr_dummy 
     2498 
     2499                      fast_wdr_dummy = fast_reservoir(ig,ib) 
     2500                      fast_reservoir(ig,ib) = MAX( (un - avail_reserve(2) * fraction_aeirrig_sw(ig) ) * fast_reservoir(ig,ib) , & 
     2501                                              fast_reservoir(ig,ib) + MIN(zero, avail_reserve(1) * fraction_aeirrig_sw(ig) * stream_reservoir(ig,ib)  - & 
     2502                                              pcent_vol_irrig * irrig_actual(ig,ib) ) ) 
     2503                      fast_wdr_dummy = fast_wdr_dummy - fast_reservoir(ig,ib) 
     2504                      irrig_fast_source(ig) = irrig_fast_source(ig) + fast_wdr_dummy 
     2505 
     2506                      stre_wdr_dummy = stream_reservoir(ig,ib) 
     2507                      stream_reservoir(ig,ib) = MAX((un - avail_reserve(1)* fraction_aeirrig_sw(ig) )*stream_reservoir(ig,ib), & 
     2508                                                stream_reservoir(ig,ib)  - & 
     2509                                                pcent_vol_irrig * irrig_actual(ig,ib) ) 
     2510                      stre_wdr_dummy = stre_wdr_dummy - stream_reservoir(ig,ib) 
     2511                      irrig_str_source(ig) = irrig_str_source(ig) + stre_wdr_dummy 
     2512 
     2513                      irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib) 
     2514 
     2515                      !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal 
     2516                      !We assume that the pot. requirement = Needs * fraction of area equipped for SW/GW 
     2517                      !In the case of surface. we also sustract the withdrawal from Fast/Stream, because both are 
     2518                      !  considered as surface water 
     2519                      IsFail_slow = ( ( irrig_needs(ig,ib)*(un - fraction_aeirrig_sw(ig)) ) > pot_slow_wdr_dummy ) 
     2520                      IsFail_fast = ( ( irrig_needs(ig,ib)*fraction_aeirrig_sw(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy ) 
     2521                      IsFail_stre = ( ( irrig_needs(ig,ib)*fraction_aeirrig_sw(ig) - fast_wdr_dummy ) > pot_stre_wdr_dummy ) 
     2522 
     2523                      IF( IsFail_stre ) THEN 
     2524                        Count_failure_stre(ig) = un 
     2525                      ENDIF 
     2526                      IF( IsFail_fast ) THEN 
     2527                        Count_failure_fast(ig) = un 
     2528                      ENDIF 
     2529                      IF( IsFail_slow ) THEN 
     2530                        Count_failure_slow(ig) = un 
     2531                      ENDIF 
     2532 
     2533                    ELSE IF (.NOT. old_irrig_scheme .AND. .NOT. select_source_irrig) THEN 
     2534                        ! For   irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir 
     2535                        ! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir 
     2536                        irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib) 
     2537 
     2538                        pot_slow_wdr_dummy = avail_reserve(3)*slow_reservoir(ig,ib) 
     2539                        pot_fast_wdr_dummy = avail_reserve(2)*fast_reservoir(ig,ib) 
     2540                        pot_stre_wdr_dummy = avail_reserve(1)*stream_reservoir(ig,ib) 
     2541                        IsFail_slow = .FALSE. ! 
     2542                        IsFail_fast = .FALSE. ! 
     2543                        IsFail_stre = .FALSE. ! 
     2544 
     2545                        irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),& 
     2546                             & pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy ) 
     2547 
     2548                        !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the 
     2549                        !Volume used for irrigation that comes from GW 
     2550                        ! Idem for irrig_fast_source and irrig_str_source 
     2551                        slow_wdr_dummy = slow_reservoir(ig,ib) 
     2552                        slow_reservoir(ig,ib) = MAX( (1-avail_reserve(3) )*slow_reservoir(ig,ib), slow_reservoir(ig,ib) + & 
     2553                             & MIN(zero, avail_reserve(2)*fast_reservoir(ig,ib) + MIN(zero, avail_reserve(1)*stream_reservoir(ig,ib)-irrig_actual(ig,ib)))) 
     2554                        slow_wdr_dummy = slow_wdr_dummy - slow_reservoir(ig,ib) 
     2555                        irrig_gw_source(ig) = irrig_gw_source(ig) + slow_wdr_dummy 
     2556 
     2557                        fast_wdr_dummy = fast_reservoir(ig,ib) 
     2558                        fast_reservoir(ig,ib) = MAX(  (1-avail_reserve(2) )*fast_reservoir(ig,ib) , & 
     2559                             fast_reservoir(ig,ib) + MIN(zero, avail_reserve(1)*stream_reservoir(ig,ib)-irrig_actual(ig,ib))) 
     2560                        fast_wdr_dummy = fast_wdr_dummy - fast_reservoir(ig,ib) 
     2561                        irrig_fast_source(ig) = irrig_fast_source(ig) + fast_wdr_dummy 
     2562 
     2563                        stre_wdr_dummy = stream_reservoir(ig,ib) 
     2564                        stream_reservoir(ig,ib) = MAX( (1-avail_reserve(1) )*stream_reservoir(ig,ib), stream_reservoir(ig,ib)-irrig_actual(ig,ib) ) 
     2565                        stre_wdr_dummy = stre_wdr_dummy - stream_reservoir(ig,ib) 
     2566                        irrig_str_source(ig) = irrig_str_source(ig) + stre_wdr_dummy 
     2567 
     2568                        irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib) 
     2569                        !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal 
     2570                        ! Because it follows the old scheme, we do not separate between surface/gw, but consider that 
     2571                        ! priority is given in this order: River, Fast and Slow reservoir. 
     2572                        IsFail_slow = ( ( irrig_needs(ig,ib) - stre_wdr_dummy - fast_wdr_dummy  ) > pot_slow_wdr_dummy ) 
     2573                        IsFail_fast = ( ( irrig_needs(ig,ib) - stre_wdr_dummy ) > pot_fast_wdr_dummy ) 
     2574                        IsFail_stre = ( irrig_needs(ig,ib) > pot_stre_wdr_dummy ) 
     2575 
     2576                        IF( IsFail_stre ) THEN 
     2577                          Count_failure_stre(ig) = un 
     2578                        ENDIF 
     2579                        IF( IsFail_fast ) THEN 
     2580                          Count_failure_fast(ig) = un 
     2581                        ENDIF 
     2582                        IF( IsFail_slow ) THEN 
     2583                          Count_failure_slow(ig) = un 
     2584                        ENDIF 
     2585 
     2586                    ELSE !Old irrigation scheme as in tag 2.0 
     2587                      !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the 
     2588                      !Volume used for irrigation that comes from GW 
     2589                      ! Idem for irrig_fast_source and irrig_str_source 
     2590                        irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib) 
     2591 
     2592                        pot_slow_wdr_dummy = slow_reservoir(ig,ib) 
     2593                        pot_fast_wdr_dummy = fast_reservoir(ig,ib) 
     2594                        pot_stre_wdr_dummy = stream_reservoir(ig,ib) 
     2595                        IsFail_slow = .FALSE. ! 
     2596                        IsFail_fast = .FALSE. ! 
     2597                        IsFail_stre = .FALSE. ! 
     2598                        irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),& 
     2599                             &   stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) ) 
     2600 
     2601                        slow_wdr_dummy = slow_reservoir(ig,ib) 
     2602                        slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + & 
     2603                             & MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib)))) 
     2604                        slow_wdr_dummy = slow_wdr_dummy - slow_reservoir(ig,ib) 
     2605                        irrig_gw_source(ig) = irrig_gw_source(ig) + slow_wdr_dummy 
     2606 
     2607                        fast_wdr_dummy = fast_reservoir(ig,ib) 
     2608                        fast_reservoir(ig,ib) = MAX( zero, & 
     2609                             &  fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib))) 
     2610                        fast_wdr_dummy = fast_wdr_dummy - fast_reservoir(ig,ib) 
     2611                        irrig_fast_source(ig) = irrig_fast_source(ig) + fast_wdr_dummy 
     2612 
     2613                        stre_wdr_dummy = stream_reservoir(ig,ib) 
     2614                        stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib) ) 
     2615                        stre_wdr_dummy = stre_wdr_dummy - stream_reservoir(ig,ib) 
     2616                        irrig_str_source(ig) = irrig_str_source(ig) + stre_wdr_dummy 
     2617 
     2618                        irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib) 
     2619                        !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal 
     2620                        ! Because it follows the old scheme, we do not separate between surface/gw, but consider that 
     2621                        ! priority is given in this order: River, Fast and Slow reservoir. 
     2622                        IsFail_slow = ( ( irrig_needs(ig,ib) - stre_wdr_dummy - fast_wdr_dummy  ) > pot_slow_wdr_dummy ) 
     2623                        IsFail_fast = ( ( irrig_needs(ig,ib) - stre_wdr_dummy ) > pot_fast_wdr_dummy ) 
     2624                        IsFail_stre = ( irrig_needs(ig,ib) > pot_stre_wdr_dummy ) 
     2625 
     2626                        IF( IsFail_stre ) THEN 
     2627                          Count_failure_stre(ig) = un 
     2628                        ENDIF 
     2629                        IF( IsFail_fast ) THEN 
     2630                          Count_failure_fast(ig) = un 
     2631                        ENDIF 
     2632                        IF( IsFail_slow ) THEN 
     2633                          Count_failure_slow(ig) = un 
     2634                        ENDIF 
     2635 
     2636                    ENDIF 
     2637 
     2638                ENDIF 
     2639            ENDDO 
    23362640          ! 
    23372641          ! Check if we cannot find the missing water in another basin of the same grid (stream reservoir only). 
    … …  
    23422646!> in another basin of the same grid. If there is water in the stream reservoir of this subbasin, we create some adduction 
    23432647!> from that subbasin to the one where we need it for irrigation. 
    2344 !>  
    2345           DO ib=1,nbasmax 
    2346  
    2347              stream_tot = SUM(stream_reservoir(ig,:)) 
    2348  
    2349              DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba) 
    2350                  
    2351                 fi = MAXLOC(stream_reservoir(ig,:)) 
    2352                 ib2 = fi(1) 
    2353  
    2354                 irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), stream_reservoir(ig,ib2)) 
    2355                 stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib) 
    2356                 irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib) 
    2357               
    2358                 stream_tot = SUM(stream_reservoir(ig,:)) 
    2359                  
    2360              ENDDO 
    2361               
    2362           ENDDO 
    2363           ! 
    2364        ENDDO 
    2365        ! 
    2366        ! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams 
    2367        ! which can feed irrigation 
    2368 ! 
    2369 !> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes 
    2370 !> to the one where we need it for irrigation. 
    2371        ! 
     2648!> 
     2649            DO ib=1,nbasmax 
     2650 
     2651               stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:)) 
     2652 
     2653                 DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba) 
     2654 
     2655                      fi = MAXLOC(stream_reservoir(ig,:)) 
     2656                      ib2 = fi(1) 
     2657 
     2658                      irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), a_stream_adduction * stream_reservoir(ig,ib2)) 
     2659                      stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib) 
     2660                      irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib) 
     2661 
     2662                      stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:)) 
     2663 
     2664                 ENDDO 
     2665 
     2666            ENDDO 
     2667          ! 
     2668      ENDDO 
     2669      ! 
     2670      ! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams 
     2671      ! which can feed irrigation 
     2672      ! 
     2673      !> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes 
     2674      !> to the one where we need it for irrigation. 
     2675      ! 
    23722676       IF (is_root_prc) THEN 
    23732677          ALLOCATE(irrig_deficit_glo(nbp_glo, nbasmax), stream_reservoir_glo(nbp_glo, nbasmax), & 
    … …  
    23982702                         ig2 = neighbours_g(ig,in) 
    23992703                         IF (ig2 .GT. 0 ) THEN 
    2400                             streams_around(in,:) = stream_reservoir_glo(ig2,:) 
     2704                            streams_around(in,:) = a_stream_adduction * stream_reservoir_glo(ig2,:) 
    24012705                            igrd(in) = ig2 
    24022706                         ENDIF 
    … …  
    24092713                         ib2=ff(2) 
    24102714                         ! 
    2411                          IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. stream_reservoir_glo(ig2,ib2) > zero ) THEN 
    2412                             adduction = MIN(irrig_deficit_glo(ig,ib), stream_reservoir_glo(ig2,ib2)) 
     2715                         IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. a_stream_adduction * stream_reservoir_glo(ig2,ib2) > zero ) THEN 
     2716                            adduction = MIN(irrig_deficit_glo(ig,ib), a_stream_adduction * stream_reservoir_glo(ig2,ib2)) 
    24132717                            stream_reservoir_glo(ig2,ib2) = stream_reservoir_glo(ig2,ib2) - adduction 
    24142718                            irrig_deficit_glo(ig,ib) = irrig_deficit_glo(ig,ib) - adduction 
    … …  
    24762780    pond_diag(:) =  zero 
    24772781    irrigation(:) = zero 
     2782    irrigdeficit(:) = zero ! 
     2783    irrigadduct(:) = zero ! 
    24782784    ! 
    24792785    ! 
    … …  
    24822788       DO ig=1,nbpt 
    24832789          IF (hydrodiag(ig,ib) > 0 ) THEN 
    2484              hydrographs(ig) = hydrographs(ig) + fast_flow(ig,ib) + slow_flow(ig,ib) + &  
    2485                   &  stream_flow(ig,ib)  
     2790             hydrographs(ig) = hydrographs(ig) + fast_flow(ig,ib) + slow_flow(ig,ib) + & 
     2791                  &  stream_flow(ig,ib) 
    24862792             slowflow_diag(ig) = slowflow_diag(ig) + slow_flow(ig,ib) 
    24872793          ENDIF 
    … …  
    24912797          flood_diag(ig) = flood_diag(ig) + flood_reservoir(ig,ib) 
    24922798          irrigation (ig) = irrigation (ig) + irrig_actual(ig,ib) + irrig_adduct(ig,ib) 
     2799          irrigdeficit(ig) = irrigdeficit(ig) + irrig_deficit(ig,ib)! 
     2800          irrigadduct (ig) = irrigadduct (ig) + irrig_adduct(ig,ib)! 
    24932801       ENDDO 
    24942802    ENDDO 
    … …  
    25022810       ! 
    25032811       irrigation(ig) = irrigation(ig)/totarea(ig) 
    2504        ! 
    2505        ! The three output types for the routing : endoheric basins,, rivers and  
     2812       irrigdeficit(ig) = irrigdeficit(ig)/totarea(ig)! 
     2813       irrigadduct(ig) = irrigadduct(ig)/totarea(ig)! 
     2814       irrig_gw_source(ig) =  irrig_gw_source(ig)/totarea(ig)! 
     2815       irrig_fast_source(ig) =  irrig_fast_source(ig)/totarea(ig)! 
     2816       irrig_str_source(ig) =  irrig_str_source(ig)/totarea(ig)! 
     2817 
     2818       ! 
     2819       ! The three output types for the routing : endoheric basins,, rivers and 
    25062820       ! diffuse coastal flow. 
    25072821       ! 
    … …  
    25132827    ! 
    25142828    flood_res = flood_diag + pond_diag 
    2515      
    2516  
    2517     !! Remove water from lake reservoir if it exceeds the maximum limit and distribute it  
     2829 
     2830 
     2831    !! Remove water from lake reservoir if it exceeds the maximum limit and distribute it 
    25182832    !! uniformly over all possible the coastflow gridcells 
    2519      
     2833 
    25202834    ! Calculate lake_overflow and remove it from lake_reservoir 
    25212835    DO ig=1,nbpt 
    … …  
    25412855    ! Transform from kg/grid-cell/dt_routing into m^3/grid-cell/s to match output unit of coastalflow 
    25422856    CALL xios_orchidee_send_field("lake_overflow_coast",lake_overflow_coast/mille/dt_routing) 
    2543     
     2857 
     2858    ! Counter of reservoir failure to assure irrigation 
     2859    CALL xios_orchidee_send_field("Count_failure_slow", Count_failure_slow) 
     2860    CALL xios_orchidee_send_field("Count_failure_fast", Count_failure_fast) 
     2861    CALL xios_orchidee_send_field("Count_failure_stre", Count_failure_stre) 
     2862 
    25442863 
    25452864  END SUBROUTINE routing_flow 
    … …  
    25902909       ! 
    25912910       lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig) 
    2592         
     2911 
    25932912       IF ( doswamps ) THEN 
    25942913          ! Calculate a return flow that will be extracted from the lake reservoir and reinserted in the soil in hydrol 
    … …  
    26352954    ! 
    26362955    IMPLICIT NONE 
    2637      
     2956 
    26382957!! INPUT VARIABLES 
    26392958    INTEGER(i_std), INTENT(in)      :: nbpt               !! Domain size (unitless) 
    … …  
    26522971 
    26532972!_ ================================================================================================================================ 
    2654     routing_area => routing_area_glo   
     2973    routing_area => routing_area_glo 
    26552974    topo_resid => topo_resid_glo 
    26562975    route_togrid => route_togrid_glo 
    … …  
    26602979    hydrodiag=>hydrodiag_glo 
    26612980    hydroupbasin=>hydroupbasin_glo 
    2662      
     2981 
    26632982    IF (is_root_prc) CALL routing_diagnostic(nbp_glo, index_g, lalo_g, resolution_g, contfrac_g, nbrivers_g,basinmap_g) 
    26642983 
    2665     routing_area => routing_area_loc   
     2984    routing_area => routing_area_loc 
    26662985    topo_resid => topo_resid_loc 
    26672986    route_togrid => route_togrid_loc 
    … …  
    26712990    hydrodiag=>hydrodiag_loc 
    26722991    hydroupbasin=>hydroupbasin_loc 
    2673      
     2992 
    26742993    CALL scatter(nbrivers_g,nbrivers) 
    26752994    CALL scatter(basinmap_g,basinmap) 
    26762995    CALL scatter(hydrodiag_glo,hydrodiag_loc) 
    26772996    CALL scatter(hydroupbasin_glo,hydroupbasin_loc) 
    2678         
     2997 
    26792998    CALL xios_orchidee_send_field("basinmap",basinmap) 
    26802999    CALL xios_orchidee_send_field("nbrivers",nbrivers) 
    … …  
    26923011       ENDIF 
    26933012    ENDIF 
    2694      
    2695          
     3013 
     3014 
    26963015  END SUBROUTINE routing_diagnostic_p 
    26973016 
    … …  
    26993018!! SUBROUTINE   : routing_diagnostic 
    27003019!! 
    2701 !>\BRIEF         This non-parallelized subroutine gives a diagnostic of the basins used. This produces some information  
     3020!>\BRIEF         This non-parallelized subroutine gives a diagnostic of the basins used. This produces some information 
    27023021!!               on the rivers which are being diagnosed. 
    27033022!! 
    27043023!! DESCRIPTION (definitions, functional, design, flags) : As not all rivers can be monitored in the model, we will only 
    27053024!! archive num_largest rivers. In this routine we will diagnose the num_largest largest rivers and print to the standard 
    2706 !! output the names of these basins and their area. The list of names of these largest rivers are taken from a list coded in the  
    2707 !! routine routing_names. As this standard output is not sufficient, we will also write it to a netCDF file with the routine  
     3025!! output the names of these basins and their area. The list of names of these largest rivers are taken from a list coded in the 
     3026!! routine routing_names. As this standard output is not sufficient, we will also write it to a netCDF file with the routine 
    27083027!! routing_diagncfile. It is important to keep for diagnostic the fraction of the largest basins in each grid box and keep information 
    27093028!! how they are linked one to the other. 
    … …  
    27263045    INTEGER(i_std), INTENT(in)                   :: nbpt                !! Domain size  (unitless) 
    27273046    INTEGER(i_std), INTENT(in)                   :: l_index(nbpt)       !! Indices of the points on the map (unitless) 
    2728     REAL(r_std), INTENT(in)                      :: lalo(nbpt,2)        !! Vector of latitude and longitudes (beware of the order !)  
     3047    REAL(r_std), INTENT(in)                      :: lalo(nbpt,2)        !! Vector of latitude and longitudes (beware of the order !) 
    27293048    REAL(r_std), INTENT(in)                      :: resolution(nbpt,2)  !! The size of each grid box in X and Y (m) 
    27303049    REAL(r_std), INTENT(in)                      :: contfrac(nbpt)      !! Fraction of land in each grid box (unitless;0-1) 
    … …  
    27413060    INTEGER(i_std), ALLOCATABLE, DIMENSION(:)    :: nb_pts              !! Number of points in the basin (unitless) 
    27423061    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: totarea             !! Total area of basin (m^2) 
    2743     REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: tmparea             !!  
     3062    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: tmparea             !! 
    27443063    INTEGER(i_std), ALLOCATABLE, DIMENSION(:)    :: topids              !! The IDs of the first num_largest basins (unitless) 
    27453064    CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:) :: basin_names         !! Names of the rivers (unitless) 
    … …  
    28183137    ALLOCATE(streams_resid(num_largest), stat=ier) 
    28193138    IF (ier /= 0) CALL ipslerr_p(3,'routing_diagnostic','Pb in allocate for streams_resid','','') 
    2820      
     3139 
    28213140    ALLOCATE(lbasin_area(num_largest,nbpt), lbasin_uparea(num_largest,nbpt), lrivercode(num_largest,nbpt), stat=ier) 
    28223141    IF (ier /= 0) CALL ipslerr_p(3,'routing_diagnostic','Pb in allocate for lbasin_area','','') 
    2823      
     3142 
    28243143    IF ( .NOT. is_root_prc) THEN 
    28253144       WRITE(numout,*) "routing_diagnostic is not suitable for running in parallel" 
    … …  
    28283147       CALL ipslerr_p(3,'routing_diagnostic','This routine is not suitable for running in parallel','','') 
    28293148    ENDIF 
    2830      
    2831      
     3149 
     3150 
    28323151    !Config Key   = RIVER_DESC 
    28333152    !Config Desc  = Writes out a description of the rivers 
    28343153    !Config If    = RIVER_ROUTING 
    28353154    !Config Def   = n 
    2836     !Config Help  = This flag allows to write out a file containing the list of  
     3155    !Config Help  = This flag allows to write out a file containing the list of 
    28373156    !Config         rivers which are beeing simulated. It provides location of outflow 
    28383157    !Config         drainage area, name and ID. 
    … …  
    28803199       ENDDO 
    28813200    ENDDO 
    2882      
     3201 
    28833202    nb_small = MIN(nb_small, 349) 
    2884      
     3203 
    28853204    ALLOCATE(basin_names(nb_small), stat=ier) 
    28863205    IF (ier /= 0) CALL ipslerr_p(3,'routing_diagnostic','Pb in allocate for basins_names','','') 
    … …  
    29223241             longest_river = MAX(longest_river, ic) 
    29233242             ! 
    2924              ! Now that we have an outflow check if it is one of the num_largest rivers.  
     3243             ! Now that we have an outflow check if it is one of the num_largest rivers. 
    29253244             ! In this case we keeps the location so we diagnose it. 
    29263245             ! 
    … …  
    31233442       tuparea(1:ii) = lbasin_uparea(idbas,1:ii) 
    31243443       ! 
    3125        CALL routing_diagcode(ii, tpts, tptbas, tuparea, tslen, MAXVAL(tslen), allstreams, upstreamchange, tcode)   
     3444       CALL routing_diagcode(ii, tpts, tptbas, tuparea, tslen, MAXVAL(tslen), allstreams, upstreamchange, tcode) 
    31263445       ! 
    31273446       lrivercode(idbas,:) = 0 
    … …  
    31783497          ENDIF 
    31793498       ENDDO 
    3180         
     3499 
    31813500       IF ( name_found > 1 ) THEN 
    31823501          DO ic=num_largest,1,-1 
    … …  
    31993518          ENDDO 
    32003519       ENDIF 
    3201         
     3520 
    32023521    ENDDO 
    32033522    ! 
    … …  
    32813600       WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. zero) 
    32823601    END IF 
    3283      
     3602 
    32843603    DEALLOCATE(pts) 
    32853604    DEALLOCATE(outpt) 
    … …  
    33013620       WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. 0.) 
    33023621    END IF 
    3303      
     3622 
    33043623  END SUBROUTINE routing_diagnostic 
    33053624  ! 
    … …  
    33083627!! 
    33093628!>\BRIEF       This subroutine determines the code in the Pfafstetter system for all points 
    3310 !!              within the given catchment.   
     3629!!              within the given catchment. 
    33113630!! 
    33123631!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    33223641!_ ================================================================================================================================ 
    33233642 
    3324   SUBROUTINE routing_diagcode(ip, tpts, tpbas, tuparea, tslen, ls, allstreams, upstreamchange, streamcode)   
     3643  SUBROUTINE routing_diagcode(ip, tpts, tpbas, tuparea, tslen, ls, allstreams, upstreamchange, streamcode) 
    33253644    ! 
    33263645    IMPLICIT NONE 
    … …  
    33883707             IF ( streamcode(ic) == indsubbas(ib) ) THEN 
    33893708                it =it+1 
    3390                 iw(it)=ic  
     3709                iw(it)=ic 
    33913710             ENDIF 
    33923711          ENDDO 
    … …  
    34073726             ENDDO 
    34083727          ELSE 
    3409              !  
     3728             ! 
    34103729             ! Else do the Pfafstetter numbering 
    3411              !  
     3730             ! 
    34123731             ! 
    34133732             ! Where do we have the 4 largest change in upstream area on this stream. 
    … …  
    34273746             ENDDO 
    34283747             ! 
    3429              ! Find all streams which are identical up to that junction and increase their code accordingly  
     3748             ! Find all streams which are identical up to that junction and increase their code accordingly 
    34303749             ! 
    34313750             DO i=1,it 
    … …  
    34743793!! 
    34753794!>\BRIEF         This subroutine creates a netCDF file containing all the informations 
    3476 !!                on the largest rivers which can be used for a refined analysis.  
     3795!!                on the largest rivers which can be used for a refined analysis. 
    34773796!! 
    34783797!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    35713890               &         TRIM(river_file_name),'(Solution ?)') 
    35723891       ENDIF 
    3573         
     3892 
    35743893       iret = NF90_DEF_DIM(fid, 'y', jjm_g, dims(2)) 
    35753894       IF (iret /= NF90_NOERR) THEN 
    … …  
    36243943    ! 1.3 Add attributes to the coordinate variables 
    36253944    ! 
    3626     iret = NF90_PUT_ATT(fid, nlonid, 'units', "degrees_east")  
     3945    iret = NF90_PUT_ATT(fid, nlonid, 'units', "degrees_east") 
    36273946    IF (iret /= NF90_NOERR) THEN 
    36283947       CALL ipslerr_p (3,'routing_diagncfile', 'Could not add attribut to variable '//lon_name//' for the file :', & 
    … …  
    37544073                ELSE 
    37554074                   basinfrac(i,j) = MIN(basinfrac(i,j), un) 
    3756                 ENDIF  
     4075                ENDIF 
    37574076             ENDDO 
    37584077          ENDDO 
    37594078          ! 
    3760           !  
     4079          ! 
    37614080          ! 2.2 Define the variables in the netCDF file 
    37624081          ! 
    … …  
    38344153          ENDIF 
    38354154          ! 
    3836           ! Longitude of outflow point  
     4155          ! Longitude of outflow point 
    38374156          att_str='Longitude_of_outflow_point' 
    38384157          iret = NF90_PUT_ATT(fid, varid, att_str, lalo(outpt(ib,1),2)) 
    … …  
    38774196          ENDIF 
    38784197          ! 
    3879           ! Average number of hops to go to the ocean for any stream  
     4198          ! Average number of hops to go to the ocean for any stream 
    38804199          att_str='Average_number_of_hops_to_ocean_for_any_stream' 
    38814200          iret = NF90_PUT_ATT(fid, varid, att_str, streams_avehops(ib)) 
    … …  
    38854204          ENDIF 
    38864205          ! 
    3887           ! Maximum number of hops to go to the ocean for any stream  
     4206          ! Maximum number of hops to go to the ocean for any stream 
    38884207          att_str='Maximum_number_of_hops_to_ocean_for_any_stream' 
    38894208          iret = NF90_PUT_ATT(fid, varid, att_str, streams_maxhops(ib)) 
    … …  
    40274346!! SUBROUTINE   : routing_basins_p 
    40284347!! 
    4029 !>\BRIEF        This parallelized subroutine computes the routing map if needed.  
     4348!>\BRIEF        This parallelized subroutine computes the routing map if needed. 
    40304349!! 
    40314350!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    40564375!    INTEGER(i_std)    :: neighbours_tmp(nbpt,8) 
    40574376!    INTEGER(i_std) :: i,j 
    4058      
     4377 
    40594378!    DO i=1,nbp_loc 
    40604379!      DO j=1,NbNeighb 
    … …  
    40634382!       ELSE 
    40644383!         neighbours_tmp(i,j)=neighbours(i,j)+nbp_para_begin(mpi_rank)-1 
    4065 !       ENDIF   
     4384!       ENDIF 
    40664385!      ENDDO 
    40674386!    ENDDO 
    40684387 
    4069     routing_area => routing_area_glo   
     4388    routing_area => routing_area_glo 
    40704389    topo_resid => topo_resid_glo 
    40714390    route_togrid => route_togrid_glo 
    … …  
    40734392    route_nbintobas => route_nbintobas_glo 
    40744393    global_basinid => global_basinid_glo 
    4075   
     4394 
    40764395    IF (is_root_prc) CALL routing_basins(nbp_glo,lalo_g, neighbours_g, resolution_g, contfrac_g) 
    40774396 
    4078     routing_area => routing_area_loc   
     4397    routing_area => routing_area_loc 
    40794398    topo_resid => topo_resid_loc 
    40804399    route_togrid => route_togrid_loc 
    … …  
    40894408    CALL scatter(route_nbintobas_glo,route_nbintobas_loc) 
    40904409    CALL scatter(global_basinid_glo,global_basinid_loc) 
    4091      
     4410 
    40924411  END SUBROUTINE routing_basins_p 
    4093   !  
    4094   
     4412  ! 
     4413 
    40954414!! ================================================================================================================================ 
    40964415!! SUBROUTINE   : routing_basins 
    40974416!! 
    40984417!>\BRIEF        This non-parallelized subroutine reads in the map of basins and flow direction to construct 
    4099 !!              the catchments of each grid box.  
     4418!!              the catchments of each grid box. 
    41004419!! 
    41014420!! DESCRIPTION (definitions, functional, design, flags) : 
    4102 !! The work is done in a number of steps which are performed locally on the  
     4421!! The work is done in a number of steps which are performed locally on the 
    41034422!! GCM grid: 
    41044423!!  1) First we find the grid-points of the high resolution routing grid which are 
    … …  
    41304449    INTEGER(i_std), INTENT(in)                    :: nbpt                  !! Domain size (unitless) 
    41314450    REAL(r_std), INTENT(in)                       :: lalo(nbpt,2)          !! Vector of latitude and longitudes (beware of the order !) 
    4132     INTEGER(i_std), INTENT(in)                    :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point  
     4451    INTEGER(i_std), INTENT(in)                    :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point 
    41334452                                                                           !! (1=North and then cloxkwise) 
    41344453    REAL(r_std), INTENT(in)                       :: resolution(nbpt,2)    !! The size of each grid box in X and Y (m) 
    … …  
    42124531    ! 
    42134532    IF (debug) THEN 
    4214        IF (ANY(diagbox .GT. nbpt)) THEN  
     4533       IF (ANY(diagbox .GT. nbpt)) THEN 
    42154534          WRITE(numout,*) "Debug diganostics : nbpt, diagbox", nbpt, diagbox 
    42164535          call ipslerr_p(3,'routing_basin', & 
    4217                &      'Problem with diagbox in debug mode.', &  
     4536               &      'Problem with diagbox in debug mode.', & 
    42184537               &      'diagbox values can''t be greater than land points number.', & 
    42194538               &      '(decrease diagbox wrong value)') 
    … …  
    42304549    !Config Def   = routing.nc 
    42314550    !Config Help  = The file provided here should alow the routing module to 
    4232     !Config         read the high resolution grid of basins and the flow direction  
     4551    !Config         read the high resolution grid of basins and the flow direction 
    42334552    !Config         from one mesh to the other. 
    42344553    !Config Units = [FILE] 
    … …  
    42804599    !! Basins : Provides a uniqe ID for each basin. These IDs are also used to get 
    42814600    !! the name of the basin from the table in routine routing_names. 
    4282     !!  
     4601    !! 
    42834602    !! Topoind :  is the topographic index for the retention time of the water in the 
    4284     !! grid box. It has been computed with the following formula : 1000 x sqrt(d^3/Dz)  
     4603    !! grid box. It has been computed with the following formula : 1000 x sqrt(d^3/Dz) 
    42854604    !! where d is the distance of the river from the current grid box to the next one 
    42864605    !! as indicated by the variable trip. 
    … …  
    43334652          ! Resolution in longitude 
    43344653          ! 
    4335           coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos )      
     4654          coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos ) 
    43364655          IF ( ip .EQ. 1 ) THEN 
    43374656             resol_lu(ip,jp,1) = ABS( lon_rel(ip+1,jp) - lon_rel(ip,jp) ) * pi/180. * R_Earth * coslat 
    … …  
    43614680    callsign = "routing_basins" 
    43624681    ok_interpol = .FALSE. 
    4363     !   
     4682    ! 
    43644683    nix=INT(MAXVAL(resolution_g(:,1))/MAXVAL(resol_lu(:,:,1)))+2 
    43654684    njx=INT(MAXVAL(resolution_g(:,2))/MAXVAL(resol_lu(:,:,2)))+2 
    … …  
    44544773    ALLOCATE (nbcoastal(nbpt), coastal_basin(nbpt,nwbas), stat=ALLOC_ERR) 
    44554774    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_basins','Pb in allocate for nbcoastal','','') 
    4456      
     4775 
    44574776    !    Order all sub points in each grid_box and find the sub basins 
    44584777    ! 
    … …  
    45274846            & nb_basin, basin_inbxid, basin_sz, basin_pts, basin_bxout, coast_pts, nwbas, basin_count,& 
    45284847            & basin_area, basin_hierarchy, basin_topoind, basin_id, basin_flowdir, outflow_grid,& 
    4529             & nbcoastal, coastal_basin)  
     4848            & nbcoastal, coastal_basin) 
    45304849       ! 
    45314850       ! 
    … …  
    45484867         & basin_hierarchy, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, & 
    45494868         & nbcoastal, coastal_basin, invented_basins) 
    4550     !  
     4869    ! 
    45514870    ! 
    45524871    IF (printlev>=1) WRITE(numout,*) 'The maximum number of basins in any grid :', MAXVAL(basin_count) 
    … …  
    47805099    !  We simplify the outflow. We only need the direction normal to the 
    47815100    !     box boundary and the 4 corners. 
    4782     !  
     5101    ! 
    47835102    ! Northern border 
    47845103    IF ( trip_bx(1,1) .EQ. 102 ) trip_bx(1,1) = 101 
    … …  
    48045123    DO jp=2,nbj-1 
    48055124       IF ( trip_bx(1,jp) .EQ. 106 .OR. trip_bx(1,jp) .EQ. 108 ) trip_bx(1,jp) = 107 
    4806     ENDDO        
     5125    ENDDO 
    48075126    ! 
    48085127    ! 
    … …  
    48545173       IF ( coords(ipos+1) /= undef_sechiba) THEN 
    48555174         IF ( COUNT(coords(ipos:nb_out) == coords(ipos)) > 1 ) THEN 
    4856             coords(ipos:nb_out-1) = coords(ipos+1:nb_out)  
     5175            coords(ipos:nb_out-1) = coords(ipos+1:nb_out) 
    48575176            coords(nb_out:nb_in) = undef_sechiba 
    48585177            nb_out = nb_out - 1 
    … …  
    48905209       DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1) 
    48915210          ll = MINLOC(coords_tmp(:), coords_tmp /= undef_sechiba) 
    4892           ind(ipos) = ll(1)  
     5211          ind(ipos) = ll(1) 
    48935212          coords_tmp(ll(1)) = undef_sechiba 
    48945213          ipos = ipos + 1 
    … …  
    48975216       DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1) 
    48985217          ll = MAXLOC(coords_tmp(:), coords_tmp /= undef_sechiba) 
    4899           ind(ipos) = ll(1)  
     5218          ind(ipos) = ll(1) 
    49005219          coords_tmp(ll(1)) = undef_sechiba 
    49015220          ipos = ipos + 1 
    … …  
    50665385          coast_pts(sz(ipb)) = bname(trans(ip)) 
    50675386          sz(trans(ip)) = 0 
    5068           pts(ipb, sz(ipb), 1) = pts(trans(ip), 1, 1)  
    5069           pts(ipb, sz(ipb), 2) = pts(trans(ip), 1, 2)  
     5387          pts(ipb, sz(ipb), 1) = pts(trans(ip), 1, 1) 
     5388          pts(ipb, sz(ipb), 2) = pts(trans(ip), 1, 2) 
    50705389       ENDDO 
    50715390    ENDIF 
    … …  
    51165435          ENDIF 
    51175436          ! 
    5118           ! Third issue : we have more than one outflow from the boxes. This could be  
     5437          ! Third issue : we have more than one outflow from the boxes. This could be 
    51195438          !             - flow into 2 or more neighboring GCM grids 
    51205439          !             - flow into a neighboring GCM grids and into the ocean or be a return flow (=97. =98, =99) 
    … …  
    51815500    ! 
    51825501    ! 
    5183     ! Sort the output by size of the various basins.  
     5502    ! Sort the output by size of the various basins. 
    51845503    ! 
    51855504    nb_basin = COUNT(sz(1:nbb) .GT. 0) 
    … …  
    52185537       CALL ipslerr_p(3,'routing_findbasins','Probleme less outflow points than basins','','') 
    52195538    ENDIF 
    5220      
     5539 
    52215540  END SUBROUTINE routing_findbasins 
    52225541  ! 
    … …  
    52245543!! SUBROUTINE   : routing_simplify 
    52255544!! 
    5226 !>\BRIEF         This subroutine symplifies the routing out of each basin by taking  
     5545!>\BRIEF         This subroutine symplifies the routing out of each basin by taking 
    52275546!!               out redundancies at the borders of the GCM box. 
    52285547!!               The aim is to have only one outflow point per basin and grid box. 
    5229 !!               But here we will not change the direction of the outflow.   
     5548!!               But here we will not change the direction of the outflow. 
    52305549!! 
    52315550!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    53535672                IF ( trip_flow(ip,jp,1) .EQ. outflow(ismall,1) .AND.& 
    53545673                     & trip_flow(ip,jp,2) .EQ. outflow(ismall,2) ) THEN 
    5355                    ! Now that we have found a point of the smallest sub-basin we  
     5674                   ! Now that we have found a point of the smallest sub-basin we 
    53565675                   ! look around for another sub-basin 
    53575676                   ib = 1 
    53585677                   not_found = .TRUE. 
    5359                    DO WHILE ( not_found .AND. ib .LE. itodo )  
     5678                   DO WHILE ( not_found .AND. ib .LE. itodo ) 
    53605679                      IF ( ib .NE. ill(1) ) THEN 
    53615680                         ibas = todopt(ib) 
    … …  
    54155734!! 
    54165735!>\BRIEF        This subroutine cuts the original basin which has more than one outflow 
    5417 !!              into as many subbasins as outflow directions.   
     5736!!              into as many subbasins as outflow directions. 
    54185737!! 
    54195738!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    54995818    ! 
    55005819    CALL routing_findrout(nbi, nbj, trip_tmp, basin_sz, basin_inbxid, nbout, outflow, trip_flow, outsz) 
    5501     !  
     5820    ! 
    55025821!    IF ( debug ) THEN 
    55035822!       DO ib = nb_in+1,nb 
    … …  
    55165835    !  Take out the small sub-basins. That is those which have only one grid box 
    55175836    !  This is only done if we need to save space in the number of basins. Else we 
    5518     !  can take it easy and keep diverging sub-basins for the moment.  
     5837    !  can take it easy and keep diverging sub-basins for the moment. 
    55195838    ! 
    55205839    IF ( nbbasins .GE. nbasmax ) THEN 
    … …  
    55435862                      DO ibb=1,nbout 
    55445863                         IF ( iip .EQ. outflow(ibb,1) .AND. jjp .EQ. outflow(ibb,2) ) THEN 
    5545                             outsz(ibb) = outsz(ibb)+1  
     5864                            outsz(ibb) = outsz(ibb)+1 
    55465865                            trip_flow(ip,jp,1) = outflow(ibb,1) 
    55475866                            trip_flow(ip,jp,2) = outflow(ibb,2) 
    … …  
    55845903          ENDIF 
    55855904       ENDDO 
    5586        ! A short verification  
     5905       ! A short verification 
    55875906       IF ( SUM(sz(nb_in+1:nb)) .NE. basin_sz) THEN 
    55885907          WRITE(numout,*) 'Lost some points while spliting the basin' 
    … …  
    55985917          CALL ipslerr_p(3,'routing_cutbasin','Lost some points while spliting the basin','','') 
    55995918       ENDIF 
    5600         
     5919 
    56015920       IF ( debug ) THEN 
    56025921          DO ib = nb_in+1,nb 
    … …  
    56235942!! 
    56245943!>\BRIEF        This subroutine finds, for each point, the distance to the outflow 
    5625 !!               point along the flowlines of the basin.  
     5944!!               point along the flowlines of the basin. 
    56265945!! 
    56275946!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    56916010             ENDIF 
    56926011             ! 
    5693              DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. iml*jml)  
     6012             DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. iml*jml) 
    56946013                cnt = cnt + 1 
    56956014                ntripi = ntripi + inc(trp,1) 
    … …  
    57076026                trp = NINT(trip(ntripi, ntripj)) 
    57086027             ENDDO 
    5709               
     6028 
    57106029             IF ( cnt .EQ. iml*jml) THEN 
    57116030                WRITE(numout,*) 'We could not route point (routing_findrout) :', ip, jp 
    … …  
    57176036                CALL ipslerr_p(3,'routing_hierarchy','We could not route point','','') 
    57186037             ENDIF 
    5719               
     6038 
    57206039             hierarchy(ip,jp) = topohier 
    5721               
     6040 
    57226041          ENDIF 
    57236042       ENDDO 
    … …  
    57316050!! 
    57326051!>\BRIEF        This subroutine simply computes the route to each outflow point 
    5733 !!              and returns the outflow point for each point in the basin.   
     6052!!              and returns the outflow point for each point in the basin. 
    57346053!! 
    57356054!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    57596078    INTEGER(i_std), DIMENSION(nbvmax,2), INTENT(out)        :: outflow   !! 
    57606079    INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2), INTENT(out) :: trip_flow !! 
    5761     INTEGER(i_std), INTENT(out)                             :: nbout     !!  
     6080    INTEGER(i_std), INTENT(out)                             :: nbout     !! 
    57626081    INTEGER(i_std), DIMENSION(nbvmax), INTENT(out)          :: outsz     !! 
    57636082    ! 
    … …  
    57906109    ! 
    57916110    ! 
    5792     !  Get the outflows and determine for each point to which outflow point it belong  
     6111    !  Get the outflows and determine for each point to which outflow point it belong 
    57936112    ! 
    57946113    nbout = 0 
    … …  
    58006119             outflow(nbout,1) = ip 
    58016120             outflow(nbout,2) = jp 
    5802           ENDIF  
     6121          ENDIF 
    58036122          IF ( trip(ip,jp) .GT. 0) THEN 
    58046123             trip_flow(ip,jp,1) = ip 
    … …  
    58156134             trp = trip(trip_flow(ip,jp,1), trip_flow(ip,jp,2)) 
    58166135             cnt = 0 
    5817              DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. nbi*nbj)  
     6136             DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. nbi*nbj) 
    58186137                cnt = cnt + 1 
    58196138                trip_flow(ip,jp,1) = trip_flow(ip,jp,1) + inc(trp,1) 
    … …  
    58626181!! 
    58636182!>\BRIEF        This subroutine puts the basins found for grid box in the global map. 
    5864 !!               Connection can only be made later when all information is together.  
     6183!!               Connection can only be made later when all information is together. 
    58656184!! 
    58666185!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    58886207    INTEGER(i_std), INTENT (in)                :: nbpt                   !! Domain size (unitless) 
    58896208    INTEGER(i_std), INTENT (in)                :: ib                     !! Current basin (unitless) 
    5890     INTEGER(i_std), INTENT(in)                 :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point  
     6209    INTEGER(i_std), INTENT(in)                 :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point 
    58916210                                                                         !! (1=North and then clockwise) 
    5892 !! LOCAL VARIABLES  
     6211!! LOCAL VARIABLES 
    58936212    REAL(r_std), DIMENSION(nbvmax,nbvmax)      :: area_bx                !! Area of each small box in the grid box (m^2) 
    58946213    INTEGER(i_std), DIMENSION(nbvmax,nbvmax)   :: trip_bx                !! The trip field for each of the smaller boxes (unitless) 
    … …  
    59246243       ! 
    59256244       basin_count(ib) = basin_count(ib)+1 
    5926        if (basin_count(ib) > nwbas) then  
     6245       if (basin_count(ib) > nwbas) then 
    59276246          WRITE(numout,*) 'ib=',ib 
    59286247          call ipslerr_p(3,'routing_globalize', & 
    5929                &      'Problem with basin_count : ', &  
     6248               &      'Problem with basin_count : ', & 
    59306249               &      'It is greater than number of allocated basin nwbas.', & 
    59316250               &      '(stop to count basins)') 
    5932        endif  
     6251       endif 
    59336252       basin_id(ib,basin_count(ib)) = basin_inbxid(ij) 
    59346253       ! 
    … …  
    60256344!! SUBROUTINE   : routing_linkup 
    60266345!! 
    6027 !>\BRIEF         This subroutine makes the connections between the basins and ensure global coherence.  
     6346!>\BRIEF         This subroutine makes the connections between the basins and ensure global coherence. 
    60286347!! 
    60296348!! DESCRIPTION (definitions, functional, design, flags) : 
    … …  
    60786397    INTEGER(i_std)                                 :: ff(1)                 !! 
    60796398    ! 
    6080     ! ERRORS  
     6399    ! ERRORS 
    60816400    LOGICAL                                        :: error1, error2, error3, error4, error5 !! (true/false) 
    60826401    ! 
    … …  
    60846403    LOGICAL, PARAMETER                             :: check = .TRUE.       !! (true/false) 
    60856404 
    6086 !_ ================================================================================================================================     
     6405!_ ================================================================================================================================ 
    60876406    error1=.FALSE. 
    60886407    error2=.FALSE. 
    … …  
    61336452                ! 
    61346453             ENDDO 
    6135              !  
     6454             ! 
    61366455             !  If we have more than one basin we will take the one which is lowest 
    61376456             !  in the hierarchy. 
    … …  
    61466465                      bop = sbl 
    61476466                   ELSE 
    6148                       ! The same hierarchy is allowed if both grids flow in about  
     6467                      ! The same hierarchy is allowed if both grids flow in about 
    61496468                      ! the same direction : 
    61506469                      IF ( ( MOD(basin_flowdir(inp,sbl)+1-1, 8)+1 .EQ. basin_flowdir(sp,sb)).OR. & 
    … …  
    61796498          ! 
    61806499          ! In case the outflow point was ocean or we did not find the correct basin we start to look 
    6181           ! around. We find two options for the outflow direction (dp1 & dm1) and the corresponding  
     6500          ! around. We find two options for the outflow direction (dp1 & dm1) and the corresponding 
    61826501          ! basin index (bp1 & bm1). 
    61836502          ! 
    … …  
    62026521                         bp1 = sbl 
    62036522                      ELSE 
    6204                          ! The same hierarchy is allowed if both grids flow in about  
     6523                         ! The same hierarchy is allowed if both grids flow in about 
    62056524                         ! the same direction : 
    62066525                         angle=MODULO(basin_flowdir(dp1,sbl)-basin_flowdir(sp,sb)+8,8) 
    … …  
    62236542                      IF ( basin_hierarchy(sp,sb) .GT. basin_hierarchy(dm1,sbl) ) THEN 
    62246543                         bm1 = sbl 
    6225                       ELSE                          
    6226                          ! The same hierarchy is allowed if both grids flow in about  
     6544                      ELSE 
     6545                         ! The same hierarchy is allowed if both grids flow in about 
    62276546                         ! the same direction : 
    62286547                         angle=MODULO(basin_flowdir(dm1,sbl)-basin_flowdir(sp,sb)+8,8) 
    … …  
    62406559             ! First deal with the case on land. 
    62416560             ! 
    6242              ! For that we need to check if the water will be able to flow out of the grid dp1 or dm1  
     6561             ! For that we need to check if the water will be able to flow out of the grid dp1 or dm1 
    62436562             ! and not return to our current grid. If it is the current grid 
    62446563             ! then we can not do anything with that neighbour. Thus we set the 
    … …  
    62506569                IF (basin_flowdir(dp1,bp1) .GT. 0) THEN 
    62516570                   outdp1 = neighbours(dp1,basin_flowdir(dp1,bp1)) 
    6252                    IF ( outdp1 .EQ. sp ) outdp1 = undef_int  
     6571                   IF ( outdp1 .EQ. sp ) outdp1 = undef_int 
    62536572                ELSE 
    62546573                   outdp1 = nbpt + 1 
    … …  
    62706589             bop = undef_int 
    62716590             IF ( outdp1 .LT. undef_int .AND. outdm1 .LT. undef_int) THEN 
    6272                 !  
     6591                ! 
    62736592                ! In this case we let the current basin flow into the smaller one 
    62746593                ! 
    … …  
    62926611             ELSE 
    62936612                ! 
    6294                 ! Now we are at the point where none of the neighboring points is suitable  
     6613                ! Now we are at the point where none of the neighboring points is suitable 
    62956614                ! or we have a coastal point. 
    62966615                ! 
    … …  
    63016620                ELSE 
    63026621                   ! 
    6303                    ! If we are on a land point with only land neighbors but no one suitable to let the  
     6622                   ! If we are on a land point with only land neighbors but no one suitable to let the 
    63046623                   ! water flow into we have to look for a solution in the current grid box. 
    6305                    !  
     6624                   ! 
    63066625                   ! 
    63076626                   IF ( bp1 .LT. 0 .AND. bm1 .LT. 0 ) THEN 
    … …  
    64026721               & .AND. basin_flowdir(sp,sb) .GT. 0) THEN 
    64036722             ! 
    6404               
     6723 
    64056724             IF (check) & 
    64066725                  WRITE(numout,*) 'There is no reason to here, this part of the code should be dead :', sp,sb 
    … …  
    64516770    IF (error2) THEN 
    64526771       CALL ipslerr_p(3,'routing_linkup', & 
    6453             &      'In the routine which make connections between the basins and ensure global coherence,', &  
     6772            &      'In the routine which make connections between the basins and ensure global coherence,', & 
    64546773            &      'there is a problem with outflow linkup without any valid direction. Try with check=.TRUE.', & 
    64556774            &      '(Perhaps there is a problem with the grid.)') 
    … …  
    64616780    ENDIF 
    64626781    IF (error4) THEN 
    6463        WRITE(numout,*) " routing_linkup : (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) ", &  
    6464             &  " .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))",sbl,sb,basin_id(sp,sbl),bid, &  
     6782       WRITE(numout,*) " routing_linkup : (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) ", & 
     6783            &  " .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))",sbl,sb,basin_id(sp,sbl),bid, & 
    64656784            &  basin_hierarchy(sp,sb),basin_hierarchy(sp,sbl) 
    64666785       CALL ipslerr_p(3,'routing_linkup', & 
    … …  
    64716790       WRITE(numout,*) 'We could not find the basin into which we need to flow' 
    64726791       WRITE(numout,*) 'Grid point ', sp, ' and basin ', sb 
    6473        WRITE(numout,*) 'Explored neighbours :', dm1, dp1  
     6792       WRITE(numout,*) 'Explored neighbours :', dm1, dp1 
    64746793       WRITE(numout,*) 'Outflow direction :', basin_flowdir(sp,sb) 
    64756794       WRITE(numout,*) 'Outlfow grid :', outflow_grid(sp,sb) 
    … …  
    65076826!!               will know how much area is upstream. It will help decide how to procede with the 
    65086827!!               the truncation later and allow to set correctly in outflow_grid the distinction 
    6509 !!               between coastal and river flow.  
     6828!!               between coastal and river flow. 
    65106829!! 
    65116830!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    66136932    WRITE(numout,*) 'The largest FETCH :', MAXVAL(fetch_basin) 
    66146933    ! 
    6615     ! Now we set for the 'num_largest' largest basins the outflow condition as stream flow  
    6616     ! (i.e. outflow_grid = -1) and all other outflow basins are set to coastal flow  
     6934    ! Now we set for the 'num_largest' largest basins the outflow condition as stream flow 
     6935    ! (i.e. outflow_grid = -1) and all other outflow basins are set to coastal flow 
    66176936    ! (i.e. outflow_grid = -2). The return flow is not touched. 
    66186937    ! 
    … …  
    66546973!>\BRIEF         This subroutine reduces the number of basins per grid to the value chosen by the user. 
    66556974!!               It also computes the final field which will be used to route the water at the 
    6656 !!               requested truncation.   
    6657 !! 
    6658 !! DESCRIPTION (definitions, functional, design, flags) :  
     6975!!               requested truncation. 
     6976!! 
     6977!! DESCRIPTION (definitions, functional, design, flags) : 
    66596978!! Truncate if needed and find the path closest to the high resolution data. 
    66606979!! 
    66616980!! The algorithm : 
    6662 !!  
     6981!! 
    66636982!! We only go through this procedure only as many times as there are basins to take out at most. 
    66646983!! This is important as it allows the simplifications to spread from one grid to the other. 
    6665 !! The for each step of the iteration and at each grid point we check the following options for  
     6984!! The for each step of the iteration and at each grid point we check the following options for 
    66666985!! simplifying the pathways of water : 
    66676986!! 1) If the basin of a grid flows into another basin of the same grid. Kill the one which only 
    66686987!!    served as a transition 
    6669 !! 2) If in one grid box we have a number of basins which flow into the ocean as coastal flow.  
    6670 !!    We kill the smallest one and put into the largest basin. There is no need to manage many  
    6671 !!    basins going into the ocean as coastal flows.  
    6672 !! 3) If we have streams run in parallel from one gird box to the others (that is these are  
     6988!! 2) If in one grid box we have a number of basins which flow into the ocean as coastal flow. 
     6989!!    We kill the smallest one and put into the largest basin. There is no need to manage many 
     6990!!    basins going into the ocean as coastal flows. 
     6991!! 3) If we have streams run in parallel from one gird box to the others (that is these are 
    66736992!!    different basins) we will put the smaller one in the larger one. This may hapen at any 
    66746993!!    level of the flow but in theory it should propagate downstream. 
    66756994!! 4) If we have two basins with the same ID but flow into different grid boxes we sacrifice 
    6676 !!    the smallest one and route it through the largest.  
     6995!!    the smallest one and route it through the largest. 
    66776996!! 
    66786997!! Obviously if any of the options find something then we skip the rest and take out the basin.:\n 
    … …  
    67587077       ENDDO 
    67597078       ! 
    6760        ! Go through the basins which need to be truncated.        
     7079       ! Go through the basins which need to be truncated. 
    67617080       ! 
    67627081       DO ibt=1,nbtruncate 
    … …  
    67827101          ENDDO 
    67837102          ! 
    6784           ! 2) Merge two basins which flow into the ocean as coastal or return flow  
    6785           ! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and  
     7103          ! 2) Merge two basins which flow into the ocean as coastal or return flow 
     7104          ! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and 
    67867105          ! have not found anything yet! 
    67877106          ! 
    … …  
    68187137          ENDIF 
    68197138          ! 
    6820           !   3) If we have basins flowing into the same grid but different basins then we put them  
     7139          !   3) If we have basins flowing into the same grid but different basins then we put them 
    68217140          !   together. Obviously we first work with the grid which has most streams running into it 
    68227141          !   and putting the smallest in the largests catchments. 
    … …  
    69477266       ENDDO 
    69487267       ! 
    6949        !       
     7268       ! 
    69507269    ENDDO 
    69517270    ! 
    … …  
    71037422                WRITE(numout,*) 'We should not be here as the basin flows into the pampa' 
    71047423                WRITE(numout,*) 'Last correct point :', pold, bold 
    7105                 WRITE(numout,*) 'It pointed to in the new variables :', route_togrid(pold, bold),route_tobasin(pold, bold)  
    7106                 WRITE(numout,*) 'The old variables gave :', outflow_grid(pold, bold), outflow_basin(pold, bold)  
     7424                WRITE(numout,*) 'It pointed to in the new variables :', route_togrid(pold, bold),route_tobasin(pold, bold) 
     7425                WRITE(numout,*) 'The old variables gave :', outflow_grid(pold, bold), outflow_basin(pold, bold) 
    71077426                WRITE(numout,*) 'Where we ended up :', igrif,ibasf 
    71087427                CALL ipslerr_p(3,'routing_truncate','Problem with routing..','','') 
    … …  
    71397458    DO ib=1, pickmax 
    71407459       ff = MAXLOC(floflo) 
    7141        ! tdo - To take into account rivers that do not flow to the oceans  
     7460       ! tdo - To take into account rivers that do not flow to the oceans 
    71427461       IF ( route_tobasin(ff(1), ff(2)) .GT. nbasmax ) THEN 
    71437462!       IF ( route_tobasin(ff(1), ff(2)) .EQ. nbasmax + 2) THEN 
    … …  
    71697488!! 
    71707489!>\BRIEF        The aim of this subroutine is to kill a basin (that is put into another larger one). 
    7171 !!              When we do this we need to be careful and change all associated variables.   
     7490!!              When we do this we need to be careful and change all associated variables. 
    71727491!! 
    71737492!! DESCRIPTION (definitions, functional, design, flags) : None 
    … …  
    72327551    inf = 0 
    72337552    DO WHILE (igrif .GT. 0) 
    7234        fetch_basin(igrif,ibasf) =  fetch_basin(igrif,ibasf) + fetch_basin(ib, tokill)  
     7553       fetch_basin(igrif,ibasf) =  fetch_basin(igrif,ibasf) + fetch_basin(ib, tokill) 
    72357554       it = outflow_grid(igrif, ibasf) 
    72367555       ibasf = outflow_basin(igrif, ibasf) 
    … …  
    72497568       ibasf = outflow_basin(igrif, ibasf) 
    72507569       igrif = it 
    7251     ENDDO    
     7570    ENDDO 
    72527571    ! 
    72537572    !  Redirect the flows which went into the basin to be killed before we change everything 
    … …  
    72917610       ENDIF 
    72927611       inflow_number(ing, inb) = inflow_number(ing, inb) - 1 
    7293         
     7612 
    72947613    ENDIF 
    72957614    ! 
    … …  
    73407659    basin_count(ib) = basin_count(ib) - 1 
    73417660    ! 
    7342   END SUBROUTINE routing_killbas  
     7661  END SUBROUTINE routing_killbas 
    73437662  ! 
    73447663!! ================================================================================================================================ 
    … …  
    77418060  ! 
    77428061!! ================================================================================================================================ 
    7743 !! SUBROUTINE   : routing_irrigmap 
     8062!! SUBROUTINE   : routing_floodmap 
    77448063!! 
    77458064!>\BRIEF         This  subroutine interpolates the 0.5x0.5 degree based map of irrigated areas to the resolution of the model. 
    … …  
    77478066!! DESCRIPTION (definitions, functional, design, flags) : None 
    77488067!! 
    7749 !! RECENT CHANGE(S): None 
     8068!! RECENT CHANGE(S): Irrigated is interpolated from slowproc as irrigated_next 
    77508069!! 
    77518070!! MAIN OUTPUT VARIABLE(S): 
    … …  
    77578076!_ ================================================================================================================================ 
    77588077 
    7759 SUBROUTINE routing_irrigmap (nbpt, index, lalo, neighbours, resolution, contfrac, & 
    7760        &                       init_irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id) 
     8078SUBROUTINE routing_floodmap (nbpt, index, lalo, neighbours, resolution, contfrac, & 
     8079        &                      init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id) 
    77618080    ! 
    77628081    IMPLICIT NONE 
    … …  
    77808099    INTEGER(i_std), INTENT(in)                     :: hist_id               !! Access to history file (unitless) 
    77818100    INTEGER(i_std), INTENT(in)                     :: hist2_id              !! Access to history file 2 (unitless) 
    7782     LOGICAL, INTENT(in)                            :: init_irrig            !! Logical to initialize the irrigation (true/false) 
     8101 
    77838102    LOGICAL, INTENT(in)                            :: init_flood            !! Logical to initialize the floodplains (true/false) 
    77848103    LOGICAL, INTENT(in)                            :: init_swamp            !! Logical to initialize the swamps (true/false) 
    77858104    ! 
    77868105!! OUTPUT VARIABLES 
    7787     REAL(r_std), INTENT(out)                       :: irrigated(:)          !! Irrigated surface in each grid box (m^2) 
     8106 
    77888107    REAL(r_std), INTENT(out)                       :: floodplains(:)        !! Surface which can be inundated in each grid box (m^2) 
    77898108    REAL(r_std), INTENT(out)                       :: swamp(:)              !! Surface which can be swamp in each grid box (m^2) 
    … …  
    77918110!! LOCAL VARIABLES 
    77928111    ! Interpolation variables 
    7793     !  
     8112    ! 
    77948113    INTEGER(i_std)                                 :: nbpmax, nix, njx, fopt !! 
    77958114    CHARACTER(LEN=30)                              :: callsign              !! 
    … …  
    78078126    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)       :: latrel                !! Latitude 
    78088127    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)       :: lonrel                !! Longitude 
    7809     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)       :: irrigated_frac        !! Irrigated fraction of the grid box (unitless;0-1) 
     8128 
    78108129    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)     :: flood_fracmax         !! Maximal flooded fraction of the grid box (unitless;0-1) 
    7811     REAL(r_std)                                    :: area_irrig            !! Irrigated surface in the grid box (m^2) 
     8130 
    78128131    REAL(r_std)                                    :: area_flood(ntype)     !! Flooded surface in the grid box (m^2) 
    78138132    REAL(r_std)                                    :: resolution_1          !! temporary variable 
    … …  
    78268145    !Config         with the area in m^2 of the area irrigated within each 
    78278146    !Config         0.5 0.5 deg grid box. The map currently used is the one 
    7828     !Config         developed by the Center for Environmental Systems Research  
     8147    !Config         developed by the Center for Environmental Systems Research 
    78298148    !Config         in Kassel (1995). 
    78308149    !Config Units = [FILE] 
    … …  
    78518170    ! 
    78528171    ALLOCATE (latrel(iml,jml), STAT=ALLOC_ERR) 
    7853     IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for latrel','','') 
     8172    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for latrel','','') 
    78548173 
    78558174    ALLOCATE (lonrel(iml,jml), STAT=ALLOC_ERR) 
    7856     IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for lonrel','','') 
    7857  
    7858     ALLOCATE (irrigated_frac(iml,jml), STAT=ALLOC_ERR) 
    7859     IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for irrigated_frac','','') 
     8175    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for lonrel','','') 
    78608176 
    78618177    ALLOCATE (flood_fracmax(iml,jml,ntype), STAT=ALLOC_ERR) 
    7862     IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for flood_fracmax','','') 
     8178    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for flood_fracmax','','') 
    78638179 
    78648180    IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lonrel, latrel, lev, tml, itau, date, dt, fid) 
    … …  
    78678183    CALL bcast(latrel) 
    78688184    ! 
    7869     IF (is_root_prc) CALL flinget(fid, 'irrig', iml, jml, lml, tml, 1, 1, irrigated_frac) 
    7870     CALL bcast(irrigated_frac) 
    78718185    IF (is_root_prc) CALL flinget(fid, 'lake', iml, jml, lml, tml, 1, 1, flood_fracmax(:,:,ilake)) 
    78728186    IF (is_root_prc) CALL flinget(fid, 'dam', iml, jml, lml, tml, 1, 1, flood_fracmax(:,:,idam)) 
    … …  
    78838197    DO ip=1,iml 
    78848198       DO jp=1,jml 
    7885           ! 
    7886           IF ( irrigated_frac(ip,jp) .LT. undef_sechiba-un) THEN 
    7887              irrigated_frac(ip,jp) = irrigated_frac(ip,jp)/100. 
    7888              IF ( irrigated_frac(ip,jp) < 0.005 ) irrigated_frac(ip,jp) = zero 
    7889           ENDIF 
    7890           ! 
     8199 
    78918200          DO itype=1,ntype 
    78928201             IF ( flood_fracmax(ip,jp,itype) .LT. undef_sechiba-1.) THEN 
    … …  
    78988207       ENDDO 
    78998208    ENDDO 
    7900      
     8209 
    79018210    IF (printlev>=2) THEN 
    79028211       WRITE(numout,*) 'lonrel : ', MAXVAL(lonrel), MINVAL(lonrel) 
    79038212       WRITE(numout,*) 'latrel : ', MAXVAL(latrel), MINVAL(latrel) 
    7904        WRITE(numout,*) 'irrigated_frac : ', MINVAL(irrigated_frac, MASK=irrigated_frac .GT. 0), & 
    7905             MAXVAL(irrigated_frac, MASK=irrigated_frac .LT. undef_sechiba) 
     8213 
    79068214       WRITE(numout,*) 'flood_fracmax : ', MINVAL(flood_fracmax, MASK=flood_fracmax .GT. 0), & 
    79078215            MAXVAL(flood_fracmax, MASK=flood_fracmax .LT. undef_sechiba) 
    … …  
    79118219    ! 
    79128220    ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR) 
    7913     IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for resol_lu','','') 
     8221    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for resol_lu','','') 
    79148222 
    79158223    ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) 
    7916     IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for mask','','') 
     8224    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for mask','','') 
    79178225    mask(:,:) = 0 
    79188226 
    … …  
    79308238          ! Resolution in longitude 
    79318239          ! 
    7932           coslat = MAX( COS( latrel(ip,jp) * pi/180. ), mincos )      
     8240          coslat = MAX( COS( latrel(ip,jp) * pi/180. ), mincos ) 
    79338241          IF ( ip .EQ. 1 ) THEN 
    79348242             resol_lu(ip,jp,1) = ABS( lonrel(ip+1,jp) - lonrel(ip,jp) ) * pi/180. * R_Earth * coslat 
    … …  
    79558263    ! Some lmargin is taken. 
    79568264    ! 
    7957     callsign = 'Irrigation map' 
     8265    callsign = 'Flood map' 
    79588266    ok_interpol = .FALSE. 
    79598267    IF (is_root_prc) THEN 
    … …  
    79698277 
    79708278    ALLOCATE(irrsub_index(nbpt, nbpmax, 2), STAT=ALLOC_ERR) 
    7971     IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for irrsub_index','','') 
     8279    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for irrsub_index','','') 
    79728280    irrsub_index(:,:,:)=0 
    79738281 
    79748282    ALLOCATE(irrsub_area(nbpt, nbpmax), STAT=ALLOC_ERR) 
    7975     IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for irrsub_area','','') 
     8283    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for irrsub_area','','') 
    79768284    irrsub_area(:,:)=zero 
    79778285 
    … …  
    79828290    ! 
    79838291    WHERE (irrsub_area < 0) irrsub_area=zero 
    7984     !   
     8292    ! 
    79858293    ! Test here if not all sub_area are larger than 0 if so, then we need to increase nbpmax 
    79868294    ! 
    79878295    DO ib=1,nbpt 
    79888296       ! 
    7989        area_irrig = 0.0 
    79908297       area_flood = 0.0 
    79918298       ! 
    … …  
    79948301          ip = irrsub_index(ib, fopt, 1) 
    79958302          jp = irrsub_index(ib, fopt, 2) 
    7996           ! 
    7997           IF (irrigated_frac(ip,jp) .LT. undef_sechiba-1.) THEN 
    7998              area_irrig = area_irrig + irrsub_area(ib,fopt)*irrigated_frac(ip,jp) 
    7999           ENDIF 
    80008303          ! 
    80018304          DO itype=1,ntype 
    … …  
    80068309       ENDDO 
    80078310       ! 
    8008        ! Put the total irrigated and flooded areas in the output variables 
    8009        ! 
    8010        IF ( init_irrig ) THEN 
     8311       ! Put the total flooded areas in the output variables 
     8312       ! 
     8313       ! 
     8314       IF ( init_flood ) THEN 
    80118315          ! if we are at the poles resolution(ib,1) = 0 
    80128316          IF (resolution(ib,1) == 0) THEN 
    8013              ! use pi*resolution(ib,2) to get the disc area 
    8014              resolution_1 = pi*resolution(ib,2) 
    8015           ELSE 
    8016              resolution_1 = resolution(ib,1) 
    8017           END IF 
    8018           irrigated(ib) = MIN(area_irrig, resolution_1*resolution(ib,2)*contfrac(ib)) 
    8019           IF ( irrigated(ib) < 0 ) THEN 
    8020              WRITE(numout,*) 'We have a problem here : ', irrigated(ib)  
    8021              WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2) 
    8022              WRITE(numout,*) area_irrig 
    8023              CALL ipslerr_p(3,'routing_irrigmap','Problem with irrigated...','','') 
    8024           ENDIF 
    8025 !!$          ! Compute a diagnostic of the map. 
    8026 !!$          IF(contfrac(ib).GT.zero) THEN 
    8027 !!$             irrigmap (ib) = irrigated(ib) / ( resolution(ib,1)*resolution(ib,2)*contfrac(ib) ) 
    8028 !!$          ELSE 
    8029 !!$             irrigmap (ib) = zero 
    8030 !!$          ENDIF 
    8031           ! 
    8032        ENDIF 
    8033        ! 
    8034        IF ( init_flood ) THEN 
    8035           ! if we are at the poles resolution(ib,1) = 0 
    8036           IF (resolution(ib,1) == 0) THEN  
    80378317             ! use pi*resolution(ib,2) to get the disc area 
    80388318             resolution_1 = pi*resolution(ib,2) 
    … …  
    80438323               & resolution_1*resolution(ib,2)*contfrac(ib)) 
    80448324          IF ( floodplains(ib) < 0 ) THEN 
    8045              WRITE(numout,*) 'We have a problem here : ', floodplains(ib)  
     8325             WRITE(numout,*) 'We have a problem here : ', floodplains(ib) 
    80468326             WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2) 
    80478327             WRITE(numout,*) area_flood 
    8048              CALL ipslerr_p(3,'routing_irrigmap','Problem with floodplains..','','') 
     8328             CALL ipslerr_p(3,'routing_floodmap','Problem with floodplains..','','') 
    80498329          ENDIF 
    80508330!!$          ! Compute a diagnostic of the map. 
    … …  
    80588338       IF ( init_swamp ) THEN 
    80598339          ! if we are at the poles resolution(ib,1) = 0 
    8060           IF (resolution(ib,1) == 0) THEN  
    8061              ! use pi*resolution(ib,2) to get the disc area  
     8340          IF (resolution(ib,1) == 0) THEN 
     8341             ! use pi*resolution(ib,2) to get the disc area 
    80628342             resolution_1 = pi*resolution(ib,2) 
    80638343          ELSE 
    … …  
    80668346          swamp(ib) = MIN(area_flood(iswamp), resolution_1*resolution(ib,2)*contfrac(ib)) 
    80678347          IF ( swamp(ib) < 0 ) THEN 
    8068              WRITE(numout,*) 'We have a problem here : ', swamp(ib)  
     8348             WRITE(numout,*) 'We have a problem here : ', swamp(ib) 
    80698349             WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2) 
    80708350             WRITE(numout,*) area_flood 
    8071              CALL ipslerr_p(3,'routing_irrigmap','Problem with swamp...','','') 
     8351             CALL ipslerr_p(3,'routing_floodmap','Problem with swamp...','','') 
    80728352          ENDIF 
    80738353!!$          ! Compute a diagnostic of the map. 
    … …  
    80838363    ! 
    80848364    ! 
    8085      
     8365 
    80868366    IF (printlev>=1) THEN 
    8087        IF ( init_irrig ) WRITE(numout,*) "Diagnostics irrigated :", MINVAL(irrigated), MAXVAL(irrigated) 
    80888367       IF ( init_flood ) WRITE(numout,*) "Diagnostics floodplains :", MINVAL(floodplains), MAXVAL(floodplains) 
    80898368       IF ( init_swamp ) WRITE(numout,*) "Diagnostics swamp :", MINVAL(swamp), MAXVAL(swamp) 
    … …  
    81128391    DEALLOCATE (latrel) 
    81138392    ! 
    8114   END SUBROUTINE routing_irrigmap 
     8393  END SUBROUTINE routing_floodmap 
    81158394  ! 
    81168395!! ================================================================================================================================ 
    … …  
    81678446    REAL(r_std), SAVE          :: totw_out             !! Sum of the water flow out to the routing scheme 
    81688447!$OMP THREADPRIVATE(totw_out) 
    8169     REAL(r_std), SAVE          :: totw_return          !!  
     8448    REAL(r_std), SAVE          :: totw_return          !! 
    81708449!$OMP THREADPRIVATE(totw_return) 
    8171     REAL(r_std), SAVE          :: totw_irrig           !!  
     8450    REAL(r_std), SAVE          :: totw_irrig           !! 
    81728451!$OMP THREADPRIVATE(totw_irrig) 
    8173     REAL(r_std), SAVE          :: totw_river           !!  
     8452    REAL(r_std), SAVE          :: totw_river           !! 
    81748453!$OMP THREADPRIVATE(totw_river) 
    8175     REAL(r_std), SAVE          :: totw_coastal         !!  
     8454    REAL(r_std), SAVE          :: totw_coastal         !! 
    81768455!$OMP THREADPRIVATE(totw_coastal) 
    81778456    REAL(r_std)                :: totarea              !! Total area of basin (m^2) 
    81788457    REAL(r_std)                :: area                 !! Total area of routing (m^2) 
    8179     INTEGER(i_std)             :: ig                   !!  
     8458    INTEGER(i_std)             :: ig                   !! 
    81808459    ! 
    81818460    ! Just to make sure we do not get too large numbers ! 
    … …  
    81948473       totw_slow = zero 
    81958474       totw_lake = zero 
    8196        totw_pond = zero  
     8475       totw_pond = zero 
    81978476       totw_in = zero 
    81988477       ! 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90

    r7576 r7709  
    3535  USE routing_highres 
    3636  USE routing_simple 
     37  USE constantes_soil 
    3738 
    3839  IMPLICIT NONE 
    … …  
    100101       rest_id,     hist_id,        hist2_id,   lalo,      & 
    101102       neighbours,  resolution,     contfrac,   stempdiag, & 
     103       soiltile,    irrig_frac,     veget_max,  irrigated_next, &     
    102104       returnflow,  reinfiltration, irrigation, riverflow, & 
    103105       coastalflow, flood_frac,     flood_res ) 
     106        
    104107 
    105108 
    … …  
    118121    REAL(r_std), INTENT(in)        :: contfrac(nbpt)       !! Fraction of land in each grid box (unitless;0-1) 
    119122    REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile 
     123    REAL(r_std), INTENT(in)        :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless) 
     124    REAL(r_std), INTENT(in)        :: veget_max(nbpt,nvm)  !! Maximal fraction of vegetation (unitless;0-1) ! 
     125    REAL(r_std), INTENT(in)        :: irrigated_next (nbpt)!! Dynamic irrig. area, calculated in slowproc and passed to routing! 
     126    REAL(r_std), INTENT(in)        :: irrig_frac(nbpt)     !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE. 
     127 
    120128 
    121129    !! 0.2 Output variables 
    … …  
    153161                                 neighbours,  resolution,     contfrac,   stempdiag, & 
    154162                                 returnflow,  reinfiltration, irrigation, riverflow, & 
    155                                  coastalflow, flood_frac,     flood_res ) 
     163                                 coastalflow, flood_frac,     flood_res,  soiltile,  & 
     164                                 irrig_frac,  veget_max,      irrigated_next) 
    156165 
    157166    ELSE IF (routing_method == 'highres') THEN 
    … …  
    206215       lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, & 
    207216       drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, & 
    208        stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id) 
    209  
     217       stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id, & 
     218       soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw)  
    210219 
    211220 
    … …  
    234243    REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile 
    235244    REAL(r_std), INTENT(in)        :: reinf_slope(nbpt)    !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1) 
     245    REAL(r_std), INTENT(in)        :: root_deficit(nbpt)   !! soil water deficit 
     246    REAL(r_std), INTENT(in)        :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless) 
     247    REAL(r_std), INTENT(in)        :: irrig_frac(nbpt)     !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE. 
     248    REAL(r_std), INTENT(in)        :: irrigated_next (nbpt)!! Dynamic irrig. area, calculated in slowproc and passed to routing 
     249    REAL(r_std), INTENT(in)        :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac 
    236250 
    237251    !! 0.2 Output variables 
    … …  
    254268            lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, & 
    255269            drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, & 
    256             stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id) 
     270            stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id, & 
     271            soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw) 
    257272 
    258273    ELSE IF (routing_method=='highres') THEN 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90

    r7525 r7709  
    113113  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: roughheight    !! Effective height for roughness (m) 
    114114!$OMP THREADPRIVATE(roughheight) 
    115   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: reinf_slope    !! slope coefficient (reinfiltration)  
     115  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: reinf_slope      !! slope coefficient (reinfiltration) 
    116116!$OMP THREADPRIVATE(reinf_slope) 
     117  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: reinf_slope_soil !! slope coefficient (reinfiltration) per soil tile 
     118!$OMP THREADPRIVATE(reinf_slope_soil) 
    117119 
    118120 
    … …  
    293295  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K) 
    294296!$OMP THREADPRIVATE(temp_sol_add) 
     297  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: root_deficit !! water deficit to reach IRRIGATION SOIL MOISTURE TARGET 
     298!$OMP THREADPRIVATE(root_deficit) 
     299  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: irrig_frac   !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE. 
     300!$OMP THREADPRIVATE(irrig_frac) 
     301  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: irrigated_next         !! Dynamic irrig. area, calculated in slowproc and passed to routing 
     302                                                                            !!  @tex $(m^{-2})$ @endtex 
     303!$OMP THREADPRIVATE(irrigated_next) 
     304  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: fraction_aeirrig_sw    !! Fraction of area equipped for irrigation from surface water, of irrig_frac 
     305                                                                            !! 1.0 here corresponds to irrig_frac, not grid cell 
     306!$OMP THREADPRIVATE(fraction_aeirrig_sw) 
     307 
    295308CONTAINS 
    296309 
    … …  
    457470                              frac_nobio,    njsc,         veget_max,      fraclut,           & 
    458471                              nwdfraclut,    tot_bare_soil,totfrac_nobio,  qsintmax,          & 
    459                               temp_growth) 
     472                              temp_growth,   irrigated_next, irrig_frac,   fraction_aeirrig_sw, & 
     473                              reinf_slope_soil) 
    460474     
    461475    !! 1.4 Initialize diffusion coefficients 
    … …  
    512526            rest_id,     hist_id,        hist2_id,   lalo,      & 
    513527            neighbours,  resolution,     contfrac,   stempdiag, & 
     528            soiltile,    irrig_frac,     veget_max,  irrigated_next, & 
    514529            returnflow,  reinfiltration, irrigation, riverflow, & 
    515530            coastalflow, flood_frac,     flood_res ) 
    … …  
    747762         & tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, & 
    748763         & humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, evap_bare_lim_ns, flood_frac, flood_res, & 
    749          & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope,& 
     764         & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope_soil,& 
    750765         & rest_id, hist_id, hist2_id,& 
    751766         & contfrac, stempdiag, & 
    … …  
    754769         & grndflux,gtemp,tot_bare_soil, & 
    755770         & lambda_snow,cgrnd_snow,dgrnd_snow,frac_snow_veg,temp_sol_add, & 
    756          & mc_layh, mcl_layh, soilmoist ) 
     771         & mc_layh, mcl_layh, soilmoist, root_deficit) 
    757772 
    758773          
    … …  
    783798            & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, & 
    784799            & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, & 
    785             & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id) 
     800            & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id, & 
     801            & soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw) 
    786802    ELSE 
    787803       !! 8.2 No routing, set variables to zero 
    … …  
    807823         lai, frac_age, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, & 
    808824         rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & 
    809          co2_flux, fco2_lu, fco2_wh, fco2_ha, temp_growth, tot_bare_soil) 
     825         co2_flux, fco2_lu, fco2_wh, fco2_ha, temp_growth, tot_bare_soil, & 
     826         irrigated_next, irrig_frac, reinf_slope, reinf_slope_soil) 
    810827 
    811828 
    … …  
    947964    CALL xios_orchidee_send_field("vegetfrac",veget) 
    948965    CALL xios_orchidee_send_field("maxvegetfrac",veget_max) 
     966    CALL xios_orchidee_send_field("irrigmap_dyn",irrigated_next) 
     967    CALL xios_orchidee_send_field("aei_sw",fraction_aeirrig_sw) 
    949968    CALL xios_orchidee_send_field("nobiofrac",frac_nobio) 
    950969    CALL xios_orchidee_send_field("soiltile",soiltile) 
    … …  
    10411060    CALL xios_orchidee_send_field("vevapnu",vevapnu/dt_sechiba) 
    10421061    CALL xios_orchidee_send_field("transpir",transpir*one_day/dt_sechiba) 
     1062    CALL xios_orchidee_send_field("transpot",transpot*one_day/dt_sechiba) 
    10431063    CALL xios_orchidee_send_field("inter",vevapwet*one_day/dt_sechiba) 
    10441064    histvar(:)=zero 
    … …  
    13951415                            ks,         nvan,      avan,     mcr,        & 
    13961416                            mcs,        mcfc,      mcw,                  & 
    1397                             assim_param, frac_age) 
     1417                            assim_param, frac_age, fraction_aeirrig_sw) 
    13981418     
    13991419    IF (printlev_loc>=3) WRITE (numout,*) 'sechiba_finalize done' 
    … …  
    15531573    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope','','') 
    15541574 
     1575    ALLOCATE (reinf_slope_soil(kjpindex, nstm),stat=ier) 
     1576    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope_soil','','') ! 
     1577 
    15551578    !salma: adding soil params 
    15561579    ALLOCATE (ks(kjpindex),stat=ier) 
    … …  
    18201843    ALLOCATE (soilmoist(kjpindex, nslm),stat=ier) 
    18211844    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilmoist','','') 
     1845 
     1846 
     1847    !1.5 Irrigation related variables 
     1848    ALLOCATE (root_deficit(kjpindex),stat=ier) 
     1849    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for root_deficit','','') ! 
     1850 
     1851    ALLOCATE (irrig_frac(kjpindex),stat=ier) 
     1852    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrig_frac','','') ! 
     1853    irrigation(:) = zero 
     1854 
     1855    ALLOCATE (irrigated_next(kjpindex),stat=ier) ! 
     1856    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable irrigated_next','','') ! 
     1857 
     1858    ALLOCATE (fraction_aeirrig_sw(kjpindex),stat=ier) ! 
     1859    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fraction_aeirrig_sw','','') 
    18221860 
    18231861    !! 1.6 Initialize indexing table for the vegetation fields.  
    … …  
    19221960    IF ( ALLOCATED (njsc)) DEALLOCATE (njsc) 
    19231961    IF ( ALLOCATED (reinf_slope)) DEALLOCATE (reinf_slope) 
     1962    IF ( ALLOCATED (reinf_slope_soil)) DEALLOCATE (reinf_slope_soil) 
    19241963 
    19251964    !salma: adding soil hydraulic params 
    … …  
    20002039    IF ( ALLOCATED (mcl_layh)) DEALLOCATE (mcl_layh) 
    20012040    IF ( ALLOCATED (soilmoist)) DEALLOCATE (soilmoist) 
     2041    IF ( ALLOCATED (root_deficit)) DEALLOCATE (root_deficit) 
     2042    IF ( ALLOCATED (irrig_frac)) DEALLOCATE (irrig_frac) 
     2043    IF ( ALLOCATED (irrigated_next)) DEALLOCATE (irrigated_next) 
    20022044 
    20032045!! 3. Clear all allocated memory 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90

    r7547 r7709  
    2020!!                   as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes). 
    2121!!                   Changes in slowproc_xios_initialize, slowproc_soilt and slowproc_finalize 
     22!!                   July 2022: New irrigation scheme. Here interpolation of new maps for 
     23!!                   the irrigation scheme 
    2224!! 
    2325!! REFERENCE(S) : 
    … …  
    7476  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: veget_max_new       !! New year fraction of vegetation type (0-1, unitless) 
    7577!$OMP THREADPRIVATE(veget_max_new) 
     78  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: irrigated_new       !! New year area equipped for irrigation  (m^{-2}) 
     79!$OMP THREADPRIVATE(irrigated_new) 
    7680  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: woodharvest         !! New year wood harvest 
    7781!$OMP THREADPRIVATE(woodharvest) 
    … …  
    288292                                  frac_nobio,    njsc,         veget_max,      fraclut,           & 
    289293                                  nwdfraclut,    tot_bare_soil,totfrac_nobio,  qsintmax,          & 
    290                                   temp_growth) 
     294                                  temp_growth,   irrigated_next, irrig_frac_next, fraction_aeirrig_sw, & 
     295                                  reinf_slope_soil) 
    291296 
    292297!! 0.1 Input variables 
    … …  
    321326    REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out)    :: nwdFraclut     !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless) 
    322327    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: reinf_slope    !! slope coef for reinfiltration 
     328   REAL(r_std),DIMENSION (kjpindex, nstm), INTENT(out)     :: reinf_slope_soil  !! slope coef for reinfiltration per soil tile 
    323329    !Salma: adding soil params 
    324330    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: ks             !! Hydraulic conductivity at saturation (mm {-1}) 
    … …  
    333339    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: deadleaf_cover !! Fraction of soil covered by dead leaves (unitless) 
    334340    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)     :: qsintmax       !! Maximum water storage on vegetation from interception (mm) 
     341    REAL(r_std), DIMENSION (kjpindex), INTENT (out)        :: irrigated_next !! Dynamic irrig. area, calculated in slowproc and passed to routing 
     342    REAL(r_std), DIMENSION (kjpindex), INTENT (out)        :: irrig_frac_next!! Dynamic irrig. fraction, calculated in slowproc and passed to routing 
     343    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: fraction_aeirrig_sw  !! Fraction of area equipped for irrigation from surface water, of irrig_frac 
     344                                                                                   !! 1.0 here corresponds to fraction of irrigated area, not grid cell 
    335345 
    336346!! 0.3 Local variables 
    337     INTEGER(i_std)                                     :: jsl 
     347    INTEGER(i_std)                                     :: ji, jsl 
    338348    REAL(r_std),DIMENSION (kjpindex,nslm)              :: land_frac         !! To ouput the clay/sand/silt fractions with a vertical dim 
    339349 
    … …  
    346356         ks,  nvan, avan, mcr, mcs, mcfc, mcw, & 
    347357         veget_max, tot_bare_soil, njsc, & 
    348          height, lcanop, veget_update, veget_year) 
     358         height, lcanop, veget_update, veget_year, fraction_aeirrig_sw) 
    349359     
    350360 
    … …  
    383393    ENDIF 
    384394 
     395    !! 5.1  Dynamic irrigation maps as output to sechiba_end 
     396    irrigated_next(:) = zero 
     397    irrig_frac_next(:) = zero 
     398    IF (do_irrigation ) THEN 
     399      irrigated_next(:) = irrigated_new(:) 
     400      ! irrig_frac calculation 
     401      DO ji=1,kjpindex 
     402          IF( (resolution(ji,1)*resolution(ji,2)*contfrac(ji) ) > min_sechiba) THEN 
     403            !SUM(routing_area(ig,:)) is totarea(ig) = m2 
     404            irrig_frac_next(ji) = MIN( soiltile(ji,irrig_st) * SUM(veget_max(ji,:)) , & 
     405                irrigated_new(ji) / (resolution(ji,1)*resolution(ji,2)*contfrac(ji) ) ) 
     406 
     407            ! Soiltile(fraction of vegtot) * SUM(veget_max) [SUM(veget_max) is vegtot in hydrol, calculated 
     408            ! differently in routing] = fraction of grid cell 
     409            ! irrigated(m2)/ SUM(routing_area(ig,:) = Fraction of grid cell 
     410            ! irrig_frac is always fraction of grid cell 
     411          ENDIF 
     412      ENDDO 
     413    END IF 
     414 
     415    !! 5.2 Calculation of reinf_slope_soil to pass to hydrol 
     416    reinf_slope_soil(:,:) = zero 
     417    DO ji=1,kjpindex 
     418        reinf_slope_soil(ji,:) = reinf_slope(ji) 
     419        IF( Reinfiltr_IrrigField .AND.  irrig_frac_next(ji) > min_sechiba ) THEN 
     420          reinf_slope_soil(ji, irrig_st) = MAX(reinf_slope(ji), reinf_slope_cropParam) 
     421        ENDIF 
     422    ENDDO 
    385423 
    386424    !! 6. Output with XIOS for variables done only once per run 
    … …  
    447485       rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & 
    448486       co2_flux, fco2_lu, fco2_wh, fco2_ha, & 
    449        temp_growth, tot_bare_soil) 
     487       temp_growth, tot_bare_soil, & 
     488       irrigated_next, irrig_frac_next, reinf_slope, reinf_slope_soil) 
    450489   
    451490!! INTERFACE DESCRIPTION 
    … …  
    478517                                                                               !! Calculated in sechiba, account for vegetation cover and  
    479518                                                                               !! effective time step to obtain gpp_d    
    480      
     519    REAL(r_std),DIMENSION (kjpindex), INTENT(in)        :: reinf_slope         !! slope coef for reinfiltration 
     520 
    481521!! 0.2 Output variables  
    482522    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out)  :: co2_flux            !! CO2 flux per average ground area (gC m^{-2} dt_stomate^{-1}) 
    … …  
    486526    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: temp_growth         !! Growth temperature (°C) - Is equal to t2m_month  
    487527    REAL(r_std), DIMENSION (kjpindex), INTENT(out)      :: tot_bare_soil       !! Total evaporating bare soil fraction in the mesh 
    488      
     528    REAL(r_std), DIMENSION (kjpindex), INTENT(out)      :: irrigated_next      !! Dynamic irrig. area, calculated in slowproc and passed to routing 
     529 
    489530!! 0.3 Modified variables 
    490531    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout)     :: lai            !! Leaf area index (m^2 m^{-2}) 
    … …  
    501542    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: deadleaf_cover !! Fraction of soil covered by dead leaves (unitless) 
    502543    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout)     :: qsintmax       !! Maximum water storage on vegetation from interception (mm) 
     544    REAL(r_std), DIMENSION (kjpindex), INTENT (inout)        :: irrig_frac_next   !! Dynamic irrig. fraction, calculated in slowproc and passed to routing 
     545    REAL(r_std),DIMENSION (kjpindex, nstm), INTENT(inout)    :: reinf_slope_soil  !!  slope coef for reinfiltration per soil tile 
    503546 
    504547!! 0.4 Local variables 
    … …  
    703746     
    704747 
     748 
     749    !! 6.1. Call to interpolation of dynamic irrigation map, if time to do so (as for vegetmax interpolation) 
     750    !! Important difference: veget map is updated every veget_update years. Irrig map_pft 
     751    !! is updated every year for now 
     752    !! Put here to use updates values of soiltile and veget_max 
     753    irrigated_next(:) = zero 
     754    IF (do_irrigation .AND. irrig_map_dynamic_flag) THEN 
     755      ! Attention: veget_year already updated, but veget_update must be >0, I.E. it must read veget maps 
     756      ! Seems logic to update irrigation maps if vegetation maps are updated too 
     757      IF ( FirstTsYear) THEN 
     758        CALL slowproc_readirrigmap_dyn(kjpindex, lalo, neighbours,  resolution, contfrac,         & 
     759             irrigated_new) 
     760        DO ji=1,kjpindex 
     761           IF( (resolution(ji,1)*resolution(ji,2)*contfrac(ji) ) > min_sechiba) THEN !Multiplication is total area(ig) = m2 
     762 
     763             irrig_frac_next(ji) = MIN( soiltile(ji,irrig_st) * SUM(veget_max(ji,:)) , & 
     764                 irrigated_new(ji) / (resolution(ji,1)*resolution(ji,2)*contfrac(ji) ) ) 
     765             ! soiltile(fraction of vegtot) * SUM(veget_max)  = fraction of grid cell 
     766             ! irrigated(m2)/ grid_cell_area = Fraction of grid cell 
     767             ! irrig_frac is always fraction of grid cell 
     768           ENDIF 
     769        ENDDO 
     770      ENDIF 
     771     !! Here irrigated_next from sechiba = irrigated_new from slowproc!! 
     772     !! Dynamic irrigation maps as output to sechiba_end 
     773    ENDIF 
     774    irrigated_next(:) = irrigated_new(:) 
     775 
     776    !! 
     777    !! 6.2 Calculation of reinf_slope_soil to pass to hydrol 
     778    IF (FirstTsYear) THEN 
     779      reinf_slope_soil(:,:) = zero 
     780      DO ji=1,kjpindex 
     781          reinf_slope_soil(ji,:) = reinf_slope(ji) 
     782          IF( Reinfiltr_IrrigField .AND.  irrig_frac_next(ji) > min_sechiba ) THEN 
     783            reinf_slope_soil(ji, irrig_st) = MAX(reinf_slope(ji), reinf_slope_cropParam) 
     784          ENDIF 
     785      ENDDO 
     786    ENDIF 
     787 
    705788    !! 7. Do some basic tests on the surface fractions updated above, only if 
    706789    !!    slowproc_veget has been done (do_slow). No change of the variables.  
    … …  
    744827                                frac_nobio, veget_max, reinf_slope,          & 
    745828                                ks,  nvan, avan, mcr, mcs, mcfc, mcw,        & 
    746                                 assim_param, frac_age ) 
     829                                assim_param, frac_age, fraction_aeirrig_sw) 
    747830 
    748831!! 0.1 Input variables 
    … …  
    766849    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3}) 
    767850    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3}) 
    768     REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (in):: assim_param   !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1}) 
     851    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: fraction_aeirrig_sw    !! Fraction of area equipped for irrigation from surface water, of irrig_frac 
     852                                                                                   !! 1.0 here corresponds to fraction of irrigated area, not grid cell 
     853     REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (in):: assim_param  !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1}) 
    769854    REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT(in):: frac_age  !! Age efficacity from STOMATE for isoprene 
    770855!! 0.4 Local variables 
    … …  
    794879    CALL restput_p (rest_id, 'frac_nobio', nbp_glo, nnobio, 1, kjit, frac_nobio, 'scatter',  nbp_glo, index_g) 
    795880 
     881    IF ( do_irrigation ) THEN 
     882          CALL restput_p (rest_id, 'irrigmap_dyn', nbp_glo, 1, 1, kjit, irrigated_new, 'scatter',  nbp_glo, index_g) 
     883          IF ( select_source_irrig ) THEN 
     884                CALL restput_p (rest_id, 'fraction_aeirrig_sw', nbp_glo, 1, 1, kjit, fraction_aeirrig_sw, 'scatter',  nbp_glo, index_g) 
     885          ENDIF 
     886    ENDIF 
    796887 
    797888    DO jf = 1, nleafages 
    … …  
    867958       ks,  nvan, avan, mcr, mcs, mcfc, mcw, & 
    868959       veget_max, tot_bare_soil, njsc, & 
    869        height, lcanop, veget_update, veget_year) 
     960       height, lcanop, veget_update, veget_year, fraction_aeirrig_sw) 
    870961     
    871962    !! INTERFACE DESCRIPTION 
    … …  
    907998    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3}) 
    908999    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     1000    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: fraction_aeirrig_sw    !! Fraction of area equipped for irrigation from surface water, of irrig_frac 
     1001                                                                                     !! 1.0 here corresponds to fraction of irrig. area, not grid cell 
    9091002 
    9101003    INTEGER(i_std), DIMENSION(kjpindex), INTENT(out)      :: njsc           !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 
    … …  
    9671060    ALLOCATE(veget_max_new(kjpindex, nvm), STAT=ier) 
    9681061    IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable veget_max_new','','') 
     1062 
     1063    ! Allocation of last irrig map area in case of change 
     1064    ALLOCATE(irrigated_new(kjpindex), STAT=ier) 
     1065    IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable irrigated_new','','') 
    9691066 
    9701067    ! Allocation of wood harvest 
    … …  
    16261723    ENDIF 
    16271724     
     1725 
     1726    !! 4.3 Dynamic irrigation map 
     1727    !  If do_irrigation, it will look to the dynamical irrig. map in restart 
     1728    !  If not dynamic irrig. map, it will be set to zero. 
     1729    !  If not found in restart, it will try to interpolate the map 
     1730 
     1731    irrigated_new(:) = zero ! 
     1732    IF ( do_irrigation ) THEN 
     1733       ! It will look into restart file 
     1734       var_name = 'irrigmap_dyn' 
     1735       CALL ioconf_setatt_p('UNITS', 'm2') 
     1736       CALL ioconf_setatt_p('LONG_NAME','Dynamical area equipped for irrigation') 
     1737       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigated_new, "gather", nbp_glo, index_g) 
     1738        
     1739       ! Now, if not found in the restart, read and interpolate from file 
     1740       IF ( ALL( irrigated_new(:) .EQ. val_exp ) ) THEN 
     1741          CALL slowproc_readirrigmap_dyn(kjpindex, lalo, neighbours,  resolution, contfrac,         & 
     1742               irrigated_new) 
     1743       ENDIF 
     1744       ! irrigated_next from sechiba (int output) = irrigated_new  in slowproc_initialize. 
     1745    ENDIF 
     1746 
     1747    ! If new priorization scheme is used, it also seek into restart, or interpolate 
     1748 
     1749    fraction_aeirrig_sw(:) = un 
     1750    IF ( do_irrigation .AND. select_source_irrig) THEN 
     1751      ! It will look into restart file 
     1752      var_name = 'fraction_aeirrig_sw' 
     1753      CALL ioconf_setatt_p('UNITS', '%') 
     1754      CALL ioconf_setatt_p('LONG_NAME','Fraction of area equipped for irrigation with surface water') 
     1755      CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., fraction_aeirrig_sw, "gather", nbp_glo, index_g) 
     1756 
     1757      ! Now, if not found in the restart, read and interpolate from file 
     1758      IF ( ALL( fraction_aeirrig_sw(:) .EQ. val_exp ) ) THEN 
     1759         CALL slowproc_read_aeisw_map(kjpindex, lalo, neighbours,  resolution, contfrac,         & 
     1760              fraction_aeirrig_sw) 
     1761      ENDIF 
     1762      ! irrigated_next from sechiba (int output) = irrigated_new in slowproc_initialize. 
     1763   ENDIF 
     1764 
     1765 
    16281766    !! 5. Some calculations always done, with and without restart files 
    16291767        
    … …  
    17071845    IF (ALLOCATED (laimap)) DEALLOCATE (laimap) 
    17081846    IF (ALLOCATED (veget_max_new)) DEALLOCATE (veget_max_new) 
     1847    IF (ALLOCATED (irrigated_new)) DEALLOCATE (irrigated_new) 
    17091848    IF (ALLOCATED (woodharvest)) DEALLOCATE (woodharvest) 
    17101849    IF (ALLOCATED (frac_nobio_new)) DEALLOCATE (frac_nobio_new) 
    … …  
    43954534  END SUBROUTINE slowproc_change_frac  
    43964535 
     4536  !! ================================================================================================================================ 
     4537  !! SUBROUTINE   : slowproc_readirrigmap_dyn 
     4538  !! 
     4539  !>\BRIEF        Function to interpolate irrigation maps 
     4540  !! 
     4541  !! DESCRIPTION  : This function interpolates the irrigation maps from original resolution to simul. resolution 
     4542  !! 
     4543  !! RECENT CHANGE(S): None 
     4544  !! 
     4545  !! MAIN OUTPUT VARIABLE(S): :: irrigmap_new 
     4546  !! 
     4547  !! REFERENCE(S) : None 
     4548  !! 
     4549  !! FLOWCHART    : None 
     4550  !! \n 
     4551  !_ ================================================================================================================================ 
     4552 
     4553  SUBROUTINE slowproc_readirrigmap_dyn(nbpt, lalo, neighbours,  resolution, contfrac,         & 
     4554       irrigmap_new) 
     4555 
     4556    USE interpweight 
     4557 
     4558    IMPLICIT NONE 
     4559 
     4560    !  0.1 INPUT 
     4561    ! 
     4562    INTEGER(i_std), INTENT(in)                             :: nbpt            !! Number of points for which the data needs 
     4563                                                                              !! to be interpolated 
     4564    REAL(r_std), DIMENSION(nbpt,2), INTENT(in)             :: lalo            !! Vector of latitude and longitudes (beware of the order !) 
     4565    INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in)   :: neighbours      !! Vector of neighbours for each grid point 
     4566                                                                              !! (1=North and then clockwise) 
     4567    REAL(r_std), DIMENSION(nbpt,2), INTENT(in)             :: resolution      !! The size in km of each grid-box in X and Y 
     4568    REAL(r_std), DIMENSION(nbpt), INTENT(in)               :: contfrac        !! Fraction of continent in the grid 
     4569    ! 
     4570    !  0.2 OUTPUT 
     4571    ! 
     4572    REAL(r_std), DIMENSION(nbpt), INTENT(out)          :: irrigmap_new       !! new irrigation map in m2 per grid cell 
     4573    ! 
     4574    !  0.3 LOCAL 
     4575    ! 
     4576    ! 
     4577    CHARACTER(LEN=80) :: filename 
     4578    INTEGER(i_std) :: ib 
     4579    ! 
     4580    ! for irrigated_new case : 
     4581 
     4582    REAL(r_std), DIMENSION(nbpt)                     :: irrigref_frac    !! irrigation fractions re-dimensioned 
     4583    REAL(r_std), DIMENSION(nbpt)                         :: airrig           !! Availability of the soilcol interpolation 
     4584    REAL(r_std)                          :: vmin, vmax       !! min/max values to use for the renormalization 
     4585    CHARACTER(LEN=80)                                    :: variablename     !! Variable to interpolate 
     4586    CHARACTER(LEN=80)                                    :: lonname, latname !! lon, lat names in input file 
     4587    REAL(r_std), DIMENSION(nvm)                          :: variabletypevals !! Values for all the types of the variable 
     4588                                                                             !!   (variabletypevals(1) = -un, not used) 
     4589    CHARACTER(LEN=50)                                    :: fractype         !! method of calculation of fraction 
     4590                                                                             !!   'XYKindTime': Input values are kinds 
     4591                                                                             !!     of something with a temporal 
     4592                                                                             !!     evolution on the dx*dy matrix' 
     4593    LOGICAL                                              :: nonegative       !! whether negative values should be removed 
     4594    CHARACTER(LEN=50)                                    :: maskingtype      !! Type of masking 
     4595                                                                             !!   'nomask': no-mask is applied 
     4596                                                                             !!   'mbelow': take values below maskvals(1) 
     4597                                                                             !!   'mabove': take values above maskvals(1) 
     4598                                                                             !!   'msumrange': take values within 2 ranges; 
     4599                                                                             !!      maskvals(2) <= SUM(vals(k)) <= maskvals(1) 
     4600                                                                             !!      maskvals(1) < SUM(vals(k)) <= maskvals(3) 
     4601                                                                             !!        (normalized by maskvals(3)) 
     4602                                                                             !!   'var': mask values are taken from a 
     4603                                                                             !!     variable inside the file (>0) 
     4604    REAL(r_std), DIMENSION(3)                            :: maskvals         !! values to use to mask (according to 
     4605                                                                             !!   `maskingtype') 
     4606    CHARACTER(LEN=250)                                   :: namemaskvar      !! name of the variable to use to mask 
     4607    CHARACTER(LEN=250)                                   :: msg 
     4608 
     4609  !_ ================================================================================================================================ 
     4610 
     4611    IF (printlev_loc >= 5) PRINT *,'  In slowproc_read irrigmap_new' 
     4612 
     4613    ! 
     4614    !Config Key   = IRRIGATION_DYN_FILE 
     4615    !Config Desc  = Name of file from which the DYNAMIC irrigation fraction map is to be read 
     4616    !Config If    = IRRIG_DYN 
     4617    !Config Def   = IRRIGmap.nc 
     4618    !Config Help  = The name of the file to be opened to read an irrigation 
     4619    !Config         map is to be given here. 
     4620    !Config Units = [FILE] 
     4621    ! 
     4622    filename = 'IRRIGmap.nc' 
     4623    CALL getin_p('IRRIGATION_DYN_FILE',filename) 
     4624    variablename = 'irrig' 
     4625 
     4626    IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_read irrigmap_new: Read and interpolate " & 
     4627         // TRIM(filename) // " for variable " // TRIM(variablename) 
     4628 
     4629    ! Assigning values to vmin, vmax 
     4630    vmin = 1. 
     4631    vmax = 1. 
     4632 
     4633    variabletypevals = -un 
     4634 
     4635    !! Variables for interpweight 
     4636    ! Type of calculation of cell fractions 
     4637    fractype = 'default' 
     4638    ! Name of the longitude and latitude in the input file 
     4639    lonname = 'lon' 
     4640    latname = 'lat' 
     4641    ! Should negative values be set to zero from input file? 
     4642    nonegative = .TRUE. 
     4643    ! Type of mask to apply to the input data (see header for more details) 
     4644    maskingtype = 'nomask' 
     4645    ! Values to use for the masking 
     4646    maskvals = (/ 0.05, 0.05 , un /) 
     4647    ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used) 
     4648    namemaskvar = '' 
     4649 
     4650    CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,        & 
     4651      contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,        & 
     4652      maskvals, namemaskvar, -1, fractype, 0., 0.,                                 & 
     4653      irrigref_frac, airrig) 
     4654 
     4655    IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_read after interpweight_2Dcont' 
     4656 
     4657    IF (printlev_loc >= 5) THEN 
     4658      WRITE(numout,*)'  slowproc_read irrigmap_new before updating loop nbpt:', nbpt 
     4659    END IF 
     4660    irrigmap_new(:) = zero 
     4661    irrigref_frac(:) = irrigref_frac(:)/100. 
     4662 
     4663    DO ib=1,nbpt 
     4664      IF (irrigref_frac(ib) < 1. .AND. irrigref_frac(ib) > 0.005 ) THEN 
     4665        irrigmap_new(ib) = irrigref_frac(ib) * resolution(ib,1)*resolution(ib,2)*contfrac(ib) 
     4666      ELSEIF (irrigref_frac(ib) > 1.) THEN 
     4667        irrigmap_new(ib) = resolution(ib,1)*resolution(ib,2)*contfrac(ib) 
     4668      !ELSE  THEN irrigated_new= zero. ALready set, put to lisibility 
     4669      ENDIF 
     4670    ENDDO 
     4671 
     4672    ! Write diagnostics 
     4673    !CALL xios_orchidee_send_field("airrig",airrig) 
     4674 
     4675    IF (printlev_loc >= 3) WRITE(numout,*) '  slowproc_read irrigmap_new ended' 
     4676 
     4677  END SUBROUTINE slowproc_readirrigmap_dyn 
     4678 
     4679  !! ================================================================================================================================ 
     4680  !! SUBROUTINE   : slowproc_read_aeisw_map 
     4681  !! 
     4682  !>\BRIEF        Function to interpolate irrigation maps 
     4683  !! 
     4684  !! DESCRIPTION  : This function interpolates the maps of areas equipped for irrigation with surface water 
     4685  !!                from original resolution to simul. resolution. This is used in the new irrigation scheme, that 
     4686  !!                restrain water availability according to environmental needs and type of equippment to irrigate. 
     4687  !! 
     4688  !! RECENT CHANGE(S): None 
     4689  !! 
     4690  !! MAIN OUTPUT VARIABLE(S): :: fraction_aeirrig_sw 
     4691  !! 
     4692  !! REFERENCE(S) : None 
     4693  !! 
     4694  !! FLOWCHART    : None 
     4695  !! \n 
     4696  !_ ================================================================================================================================ 
     4697 
     4698 
     4699  SUBROUTINE slowproc_read_aeisw_map(nbpt, lalo, neighbours,  resolution, contfrac,         & 
     4700     fraction_aeirrig_sw) 
     4701 
     4702    USE interpweight 
     4703 
     4704    IMPLICIT NONE 
     4705 
     4706    !  0.1 INPUT 
     4707    ! 
     4708    INTEGER(i_std), INTENT(in)                             :: nbpt            !! Number of points for which the data needs 
     4709                                                                              !! to be interpolated 
     4710    REAL(r_std), DIMENSION(nbpt,2), INTENT(in)             :: lalo            !! Vector of latitude and longitudes (beware of the order !) 
     4711    INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in)   :: neighbours      !! Vector of neighbours for each grid point 
     4712                                                                              !! (1=North and then clockwise) 
     4713    REAL(r_std), DIMENSION(nbpt,2), INTENT(in)             :: resolution      !! The size in km of each grid-box in X and Y 
     4714    REAL(r_std), DIMENSION(nbpt), INTENT(in)               :: contfrac        !! Fraction of continent in the grid 
     4715    ! 
     4716    !  0.2 OUTPUT 
     4717    ! 
     4718    REAL(r_std), DIMENSION(nbpt), INTENT(out)          :: fraction_aeirrig_sw       !! Fraction of area equipped for irrigation from surface water, of irrig_frac 
     4719                                                                                ! 1.0 here corresponds to fraction of irrig. area, not grid cell 
     4720    ! 
     4721    !  0.3 LOCAL 
     4722    CHARACTER(LEN=80) :: filename 
     4723    INTEGER(i_std) :: ib 
     4724    ! 
     4725    ! for irrigated_new case : 
     4726 
     4727    REAL(r_std), DIMENSION(nbpt)                     :: irrigref_frac    !! irrigation fractions re-dimensioned 
     4728    REAL(r_std), DIMENSION(nbpt)                         :: airrig           !! Availability of the soilcol interpolation 
     4729    REAL(r_std)                          :: vmin, vmax       !! min/max values to use for the renormalization 
     4730    CHARACTER(LEN=80)                                    :: variablename     !! Variable to interpolate 
     4731    CHARACTER(LEN=80)                                    :: lonname, latname !! lon, lat names in input file 
     4732    REAL(r_std), DIMENSION(nvm)                          :: variabletypevals !! Values for all the types of the variable 
     4733                                                                             !!   (variabletypevals(1) = -un, not used) 
     4734    CHARACTER(LEN=50)                                    :: fractype         !! method of calculation of fraction 
     4735                                                                             !!   'XYKindTime': Input values are kinds 
     4736                                                                             !!     of something with a temporal 
     4737                                                                             !!     evolution on the dx*dy matrix' 
     4738    LOGICAL                                              :: nonegative       !! whether negative values should be removed 
     4739    CHARACTER(LEN=50)                                    :: maskingtype      !! Type of masking 
     4740                                                                             !!   'nomask': no-mask is applied 
     4741                                                                             !!   'mbelow': take values below maskvals(1) 
     4742                                                                             !!   'mabove': take values above maskvals(1) 
     4743                                                                             !!   'msumrange': take values within 2 ranges; 
     4744                                                                             !!      maskvals(2) <= SUM(vals(k)) <= maskvals(1) 
     4745                                                                             !!      maskvals(1) < SUM(vals(k)) <= maskvals(3) 
     4746                                                                             !!        (normalized by maskvals(3)) 
     4747                                                                             !!   'var': mask values are taken from a 
     4748                                                                             !!     variable inside the file (>0) 
     4749    REAL(r_std), DIMENSION(3)                            :: maskvals         !! values to use to mask (according to 
     4750                                                                             !!   `maskingtype') 
     4751    CHARACTER(LEN=250)                                   :: namemaskvar      !! name of the variable to use to mask 
     4752    CHARACTER(LEN=250)                                   :: msg 
     4753 
     4754  !_ ================================================================================================================================ 
     4755 
     4756    IF (printlev_loc >= 5) PRINT *,'  In slowproc_read irrigmap_new' 
     4757 
     4758    ! 
     4759    !Config Key   = FRACTION_AEI_SW_FILE 
     4760    !Config Desc  = Name of file with AEI with SW 
     4761    !Config If    = SELECT_SOURCE_IRRIG 
     4762    !Config Def   = AEI_SW_pct.nc 
     4763    !Config Help  = The name of the file to be opened to read an AIE from SW 
     4764    !Config         map is to be given here. 
     4765    !Config Units = [FILE] 
     4766    ! 
     4767    filename = 'AEI_SW_pct.nc' 
     4768    CALL getin_p('FRACTION_AEI_SW_FILE',filename) 
     4769    variablename = 'aeisw_pct' 
     4770 
     4771    IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_read_aeisw_map: Read and interpolate " & 
     4772         // TRIM(filename) // " for variable " // TRIM(variablename) 
     4773 
     4774    ! Assigning values to vmin, vmax 
     4775    vmin = 1. 
     4776    vmax = 1. 
     4777 
     4778    variabletypevals = -un 
     4779 
     4780    !! Variables for interpweight 
     4781    ! Type of calculation of cell fractions 
     4782    fractype = 'default' 
     4783    ! Name of the longitude and latitude in the input file 
     4784    lonname = 'lon' 
     4785    latname = 'lat' 
     4786    ! Should negative values be set to zero from input file? 
     4787    nonegative = .TRUE. 
     4788    ! Type of mask to apply to the input data (see header for more details) 
     4789    maskingtype = 'mabove' 
     4790    ! Values to use for the masking 
     4791    maskvals = (/ 0.05, 0.05 , un /) 
     4792    ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used) 
     4793    namemaskvar = '' 
     4794 
     4795    CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,        & 
     4796      contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,        & 
     4797      maskvals, namemaskvar, -1, fractype, 0., 0.,                                 & 
     4798      irrigref_frac, airrig) 
     4799 
     4800    IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_read after interpweight_2Dcont' 
     4801 
     4802    IF (printlev_loc >= 5) THEN 
     4803      WRITE(numout,*)'  slowproc_read_aeisw_map before updating loop nbpt:', nbpt 
     4804    END IF 
     4805    fraction_aeirrig_sw(:) = irrigref_frac(:)/100. 
     4806    ! Write diagnostics 
     4807    !CALL xios_orchidee_send_field("airrig",airrig) 
     4808    DO ib=1,nbpt 
     4809 
     4810      fraction_aeirrig_sw(ib) = MIN(fraction_aeirrig_sw(ib), 0.99 ) 
     4811      fraction_aeirrig_sw(ib) = MAX(fraction_aeirrig_sw(ib), 0.01 ) 
     4812 
     4813    ENDDO 
     4814 
     4815 
     4816    IF (printlev_loc >= 3) WRITE(numout,*) 'slowproc_read_aeisw_map ended' 
     4817 
     4818  END SUBROUTINE slowproc_read_aeisw_map 
     4819 
    43974820END MODULE slowproc 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml

    r7576 r7709  
    384384    <field id="zcarbwh" name="zcarbwh" long_name="Wood harvest CO2 flux (variable in interface to LMDZ)" unit="kgC m-2 s-1" /> 
    385385    <field id="zcarbha" name="zcarbha" long_name="Crop harvesting CO2 flux (variable in interface to LMDZ)" unit="kgC m-2 s-1" /> 
     386 
     387    <!-- Irrigation scheme variables                                                                  --> 
     388    <field id="irrigmap_dyn" name="irrigmap_dyn" long_name="Dyn. map of irrigated areas" unit="m^2"/> 
     389    <field id="aei_sw" name="aei_sw" long_name="Area Equipped for irrigation with surface water" unit="m^2"/> 
     390    <field id="root_deficit" name="root_deficit" long_name="Soil Deficit " unit="kg/m^2"/> 
     391    <field id="root_mc_fc" name="root_mc_fc" long_name="FC SM, for root zone " unit="kg/m^2"/> 
     392    <field id="irrig_gw_source" name="irrig_gw_source_mean" long_name="Irrigation flux from slow reservoir" unit="mm/s"/> 
     393    <field id="irrig_fast_source" name="irrig_fast_source_mean" long_name="Irrigation flux from fast reservoir" unit="mm/s"/> 
     394    <field id="irrig_str_source" name="irrig_str_source_mean" long_name="Irrigation flux from stream reservoir" unit="mm/s"/> 
     395    <field id="vbeta3pot" name="vbeta3pot" long_name="Beta for Potential Transpiration" unit="mm/d" grid_ref="grid_nvm"/> 
     396    <field id="transpot" name="transpot" long_name="Potential transpiration" unit="mm/d" grid_ref="grid_nvm"/> 
     397    <field id="irrig_deficit" name="irrig_deficit" long_name="Difference between demand and irrigation" unit="mm/s"/> 
     398    <field id="irrig_adduct" name="irrig_adduct" long_name="Irrigation from adduction in nearby basins and grid cells" unit="mm/s"/> 
     399    <field id="Count_failure_slow" name="Count_failure_slow" long_name="Counter of failure of slow res. to fulfill irrigation water" unit="-"/> 
     400    <field id="Count_failure_fast" name="Count_failure_fast" long_name="Counter of failure of fast res. to fulfill irrigation water" unit="-"/> 
     401    <field id="Count_failure_stre" name="Count_failure_stre" long_name="Counter of failure of stream res. to fulfill irrigation water" unit="-"/> 
     402 
    386403  </field_group> 
    387404 

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml

    r7613 r7709  
    343343    <field field_ref="snowtemp_read_current" grid_ref="grid_nsnow_out"  level="12"/> 
    344344    <field field_ref="mask_snow_interp_out" grid_ref="grid_nsnow_out"  level="12"/> 
     345 
     346     <!-- Irrigation scheme  --> 
     347     <field field_ref="irrigation" name="Qirrig" level="8"/> 
     348     <field field_ref="irrig_gw_source" name="Qirrig_gw" level="8"/> 
     349     <field field_ref="irrig_gw_source" name="irrig_gw_source" level="8" unit="mm/d" > this*86400  </field> 
     350     <field field_ref="irrig_fast_source" name="Qirrig_fast" level="8"/> 
     351     <field field_ref="irrig_fast_source" name="irrig_fast_source" level="8" unit="mm/d" > this*86400  </field> 
     352     <field field_ref="irrig_str_source" name="Qirrig_str" level="8"/> 
     353     <field field_ref="irrig_str_source" name="irrig_str_source" level="8" unit="mm/d" > this*86400  </field> 
     354     <field field_ref="netirrig" name="Qirrig_req" long_name="Irrigation requirement" level="8"/> 
     355     <field field_ref="soilflx" level="8"/> 
     356     <field field_ref="irrig_deficit" name="Qirrig_deficit" level="8"/> 
     357     <field field_ref="irrig_deficit" name="irrig_deficit" level="8" unit="mm/d" > this*86400  </field> 
     358     <field field_ref="irrig_adduct" name="Qirrig_adduct" level="8"/> 
     359     <field field_ref="irrig_adduct" name="irrig_adduct" level="8" unit="mm/d" > this*86400  </field> 
     360     <field field_ref="Count_failure_slow" name="Count_failure_slow" level="8"/> 
     361     <field field_ref="Count_failure_fast" name="Count_failure_fast" level="8"/> 
     362     <field field_ref="Count_failure_stre" name="Count_failure_stre" level="8"/> 
     363     <!--field field_ref="PotSurfT" level="8"/--> 
     364     <field field_ref="transpot" grid_ref="grid_nvm_out" level="8"/> 
     365     <field field_ref="vbeta3pot" grid_ref="grid_nvm_out" level="8"/> 
     366     <field field_ref="root_deficit" level="8"/> 
     367     <field field_ref="root_mc_fc" level="8"/> 
     368     <field field_ref="irrigmap_dyn" level="8"/> 
     369     <field field_ref="aei_sw" level="8"/> 
     370 
    345371   </field_group> 
    346372  </file> 
Note: See TracChangeset for help on using the changeset viewer.