Changeset 7709 for branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90
- Timestamp:
- 2022-07-20T11:30:43+02:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90
r7476 r7709 24 24 !! We can also impose kfact_root=1 in all soil layers to cancel the effect of 25 25 !! 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. 26 28 !! 27 29 !! REFERENCE(S) : … … 366 368 !! @tex $(m^{3} m^{-3})$ @endtex 367 369 !$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) 368 376 369 377 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc_read_prev !! Soil moisture from file at previous timestep in the file … … 565 573 & humrel, vegstress, drysoil_frac, evapot, evapot_penm, evap_bare_lim, evap_bare_lim_ns, & 566 574 & 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,& 568 576 & contfrac, stempdiag, & 569 577 & temp_air, pb, u, v, tq_cdrag, swnet, pgflux, & … … 571 579 & grndflux,gtemp,tot_bare_soil, & 572 580 & 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 ) 574 582 575 583 !! 0. Variable and parameter declaration … … 607 615 REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception 608 616 REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration 609 REAL(r_std),DIMENSION (kjpindex ), INTENT (in) :: reinf_slope !! Slope coef617 REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: reinf_slope_soil !! Slope coef per soil tile 610 618 611 619 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ks !! Hydraulic conductivity at saturation (mm {-1}) … … 647 655 REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_melt !! Total melt 648 656 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 650 659 !! 0.3 Modified variables 651 660 … … 794 803 !! 3.5 computes soil hydrology ==>hydrol_soil 795 804 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, & 797 806 transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 798 807 returnflow, reinfiltration, irrigation, & … … 800 809 k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,& 801 810 stempdiag,snow,snowdz, tot_bare_soil, u, v, tq_cdrag, & 802 mc_layh, mcl_layh )811 mc_layh, mcl_layh, root_deficit, veget) 803 812 804 813 ! The update fluxes come from hydrol_vegupd 805 814 drainage(:) = drainage(:) + drain_upd(:) 806 815 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) 807 820 808 821 … … 1691 1704 ALLOCATE (mc(kjpindex,nslm,nstm),stat=ier) 1692 1705 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','','') ! 1693 1712 1694 1713 … … 2365 2384 IF (ALLOCATED (vegetmax_soil)) DEALLOCATE (vegetmax_soil) 2366 2385 IF (ALLOCATED (mc)) DEALLOCATE (mc) 2386 IF (ALLOCATED (root_mc_fc)) DEALLOCATE (root_mc_fc) 2387 IF (ALLOCATED (nslm_root)) DEALLOCATE (nslm_root) 2367 2388 IF (ALLOCATED (soilmoist)) DEALLOCATE (soilmoist) 2368 2389 IF (ALLOCATED (soilmoist_liquid)) DEALLOCATE (soilmoist_liquid) … … 2728 2749 INTEGER(i_std) :: iiref !! To identify the mc_lins where k_lin and d_lin 2729 2750 !! 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 2730 2755 2731 2756 !_ ================================================================================================================================ … … 2880 2905 !! 2 Compute the root density profile if not ok_dynroot 2881 2906 !! 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 2882 2910 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 2888 2919 DO jv = 1,nvm 2889 DO jsl = 2, nslm-1 2920 DO ji=1, kjpindex 2921 2922 2890 2923 nroot(ji,jv,jsl) = (EXP(-humcste(jv)*zz(jsl)/mille)) * & 2891 2924 & (EXP(humcste(jv)*dz(jsl)/mille/deux) - & … … 2893 2926 & (EXP(-humcste(jv)*dz(2)/mille/deux) & 2894 2927 & -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 2895 2934 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 2898 2944 nroot(ji,jv,nslm) = (EXP(humcste(jv)*dz(nslm)/mille/deux) -un) * & 2899 2945 & EXP(-humcste(jv)*zz(nslm)/mille) / & … … 2902 2948 ENDDO 2903 2949 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 2904 2962 END IF 2905 2963 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 2908 2980 !- 2909 2981 !! 3 Compute the profile for a and n … … 3592 3664 !_ hydrol_soil 3593 3665 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, & 3595 3667 & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 3596 3668 & returnflow, reinfiltration, irrigation, & … … 3598 3670 & k_litt, litterhumdiag, humrel,vegstress, drysoil_frac, & 3599 3671 & 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) 3601 3673 ! 3602 3674 ! interface description … … 3608 3680 INTEGER(i_std), INTENT(in) :: kjpindex 3609 3681 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 3610 3683 INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class 3611 3684 !! in the grid cell (1-nscm, unitless) … … 3622 3695 REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: transpir !! Transpiration 3623 3696 !! @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex 3624 REAL(r_std), DIMENSION (kjpindex ), INTENT (in) :: reinf_slope !! Fraction of surface runoff that reinfiltrates3697 REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: reinf_slope_soil !! Fraction of surface runoff that reinfiltrates per soil tile 3625 3698 !! (unitless, [0-1]) 3626 3699 REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: returnflow !! Water returning to the soil from the bottom … … 3672 3745 REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out) :: mcl_layh !! Volumetric liquid water content for each soil layer 3673 3746 !! 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 3675 3750 !! 0.3 Modified variables 3676 3751 … … 3698 3773 REAL(r_std), DIMENSION(kjpindex) :: reinfiltration_soil !! Water from the routing back to the top of the 3699 3774 !! soil @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex 3700 REAL(r_std), DIMENSION(kjpindex ) :: irrigation_soil !! Water from irrigation returning to soil moisture3775 REAL(r_std), DIMENSION(kjpindex,nstm) :: irrigation_soil !! Water from irrigation returning to soil moisture per soil tile 3701 3776 !! @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex 3702 3777 REAL(r_std), DIMENSION(kjpindex) :: flux_infilt !! Water to infiltrate … … 3766 3841 REAL(r_std), DIMENSION (kjpindex) :: check_over !! Water conservation diagnostic at routine scale 3767 3842 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 - 3769 3845 ! Diagnostic of the vertical soil water fluxes 3770 3846 REAL(r_std), DIMENSION (kjpindex,nslm) :: qflux !! Local upward flux into soil layer … … 3796 3872 returnflow_soil(:) = zero 3797 3873 reinfiltration_soil(:) = zero 3798 irrigation_soil(: ) = zero3874 irrigation_soil(:,:) = zero 3799 3875 qflux_ns(:,:,:) = zero 3800 3876 mc_layh(:,:) = zero ! for thermosoil … … 3838 3914 ! AD16*** The transformation in 0.2 and 0.3 is likely to induce conservation problems 3839 3915 ! 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 3851 3932 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 3853 3954 3854 3955 !! -- START MAIN LOOP (prognostic loop to update mc and mcl) OVER SOILTILES … … 3901 4002 ! Here, temp is used as a temporary variable to store the min of water to infiltrate vs evaporate 3902 4003 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) & 3904 4005 - MIN(ae_ns(ji,jst),zero) - MIN(subsinksoil(ji),zero) + precisol_ns(ji,jst), & 3905 4006 MAX(ae_ns(ji,jst),zero) + MAX(subsinksoil(ji),zero) ) … … 3912 4013 ! - eventually, water2infilt holds all fluxes to the soil surface except precisol (reduced by water2extract) 3913 4014 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) & 3915 4017 - MIN(ae_ns(ji,jst),zero) - MIN(subsinksoil(ji),zero) + precisol_ns(ji,jst) & 3916 4018 - temp(ji) … … 3955 4057 IF ( .NOT. doponds ) THEN ! this is the general case... 3956 4058 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) 3958 4060 ENDDO 3959 4061 ELSE … … 4341 4443 tmcs_litter(:) = tmcs_litter(:) + sms(:,jsl) 4342 4444 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 4344 4488 ! Soil wetness profiles (W-Ww)/(Ws-Ww) 4345 4489 ! soil_wet_ns is the ratio of available soil moisture to max available soil moisture
Note: See TracChangeset
for help on using the changeset viewer.
