Changeset 6954 for branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90
- Timestamp:
- 2020-12-03T18:16:43+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90
r6189 r6954 19 19 !! This routine was originaly developed by Patricia deRosnay. 20 20 !! 21 !! RECENT CHANGE(S) : None 22 !! 21 !! RECENT CHANGE(S) : November 2020: It is possible to define soil hydraulic parameters from maps, 22 !! as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes). 23 !! Here, it leads to change dimensions and indices. 24 !! We can also impose kfact_root=1 in all soil layers to cancel the effect of 25 !! roots on ks profile (keyword KFACT_ROOT_TYPE). 26 !! 23 27 !! REFERENCE(S) : 24 28 !! - de Rosnay, P., J. Polcher, M. Bruen, and K. Laval, Impact of a physically based soil … … 43 47 !! of a new soil freezing scheme for a land-surface model with physically-based hydrology. 44 48 !! The Cryosphere, 6, 407-430, doi: 10.5194/tc-6-407-2012. \n 49 !! - Tafasca S. (2020). Evaluation de lâimpact des propriétés du sol sur lâhydrologie simulee dans le 50 !! modÚle ORCHIDEE, PhD thesis, Sorbonne Universite. \n 45 51 !! 46 52 !! SVN : … … 91 97 ! one dimension array allocated, computed, saved and got in hydrol module 92 98 ! Values per soil type 93 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: nvan !! Van Genuchten coeficients n (unitless) 94 ! RK: 1/n=1-m 95 !$OMP THREADPRIVATE(nvan) 96 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: avan !! Van Genuchten coeficients a 97 !! @tex $(mm^{-1})$ @endtex 98 !$OMP THREADPRIVATE(avan) 99 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcr !! Residual volumetric water content 100 !! @tex $(m^{3} m^{-3})$ @endtex 101 !$OMP THREADPRIVATE(mcr) 102 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcs !! Saturated volumetric water content 103 !! @tex $(m^{3} m^{-3})$ @endtex 104 !$OMP THREADPRIVATE(mcs) 105 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: ks !! Hydraulic conductivity at saturation 106 !! @tex $(mm d^{-1})$ @endtex 107 !$OMP THREADPRIVATE(ks) 99 108 100 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: pcent !! Fraction of saturated volumetric soil moisture above 109 101 !! which transpir is max (0-1, unitless) 110 !$OMP THREADPRIVATE(pcent) 111 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcfc !! Volumetric water content at field capacity 112 !! @tex $(m^{3} m^{-3})$ @endtex 113 !$OMP THREADPRIVATE(mcfc) 114 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcw !! Volumetric water content at wilting point 115 !! @tex $(m^{3} m^{-3})$ @endtex 116 !$OMP THREADPRIVATE(mcw) 102 !$OMP THREADPRIVATE(pcent) 117 103 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mc_awet !! Vol. wat. cont. above which albedo is cst 118 104 !! @tex $(m^{3} m^{-3})$ @endtex … … 216 202 !$OMP THREADPRIVATE(kfact_root) 217 203 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: kfact !! Factor to reduce Ks with depth (unitless) 218 !! DIM = nslm * nscm204 !! DIM = nslm * kjpindex 219 205 !$OMP THREADPRIVATE(kfact) 220 206 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: zz !! Depth of nodes [znh in vertical_soil] transformed into (mm) … … 228 214 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_lin !! 50 Vol. Wat. Contents to linearize K and D, for each texture 229 215 !! @tex $(m^{3} m^{-3})$ @endtex 230 !! DIM = imin:imax * nscm216 !! DIM = imin:imax * kjpindex 231 217 !$OMP THREADPRIVATE(mc_lin) 232 218 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: k_lin !! 50 values of unsaturated K, for each soil layer and texture 233 219 !! @tex $(mm d^{-1})$ @endtex 234 !! DIM = imin:imax * nslm * nscm220 !! DIM = imin:imax * nslm * kjpindex 235 221 !$OMP THREADPRIVATE(k_lin) 236 222 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: d_lin !! 50 values of diffusivity D, for each soil layer and texture 237 223 !! @tex $(mm^2 d^{-1})$ @endtex 238 !! DIM = imin:imax * nslm * nscm224 !! DIM = imin:imax * nslm * kjpindex 239 225 !$OMP THREADPRIVATE(d_lin) 240 226 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: a_lin !! 50 values of the slope in K=a*mc+b, for each soil layer and texture 241 227 !! @tex $(mm d^{-1})$ @endtex 242 !! DIM = imin:imax * nslm * nscm228 !! DIM = imin:imax * nslm * kjpindex 243 229 !$OMP THREADPRIVATE(a_lin) 244 230 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: b_lin !! 50 values of y-intercept in K=a*mc+b, for each soil layer and texture 245 231 !! @tex $(m^{3} m^{-3})$ @endtex 246 !! DIM = imin:imax * nslm * nscm232 !! DIM = imin:imax * nslm * kjpindex 247 233 !$OMP THREADPRIVATE(b_lin) 248 234 … … 259 245 !$OMP THREADPRIVATE(kk) 260 246 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: avan_mod_tab !! VG parameter a modified from exponantial profile 261 !! @tex $(mm^{-1})$ @endtex !! DIMENSION (nslm, nscm)247 !! @tex $(mm^{-1})$ @endtex !! DIMENSION (nslm,kjpindex) 262 248 !$OMP THREADPRIVATE(avan_mod_tab) 263 249 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: nvan_mod_tab !! VG parameter n modified from exponantial profile 264 !! (unitless) !! DIMENSION (nslm, nscm)250 !! (unitless) !! DIMENSION (nslm,kjpindex) 265 251 !$OMP THREADPRIVATE(nvan_mod_tab) 266 252 … … 450 436 !_ ================================================================================================================================ 451 437 452 SUBROUTINE hydrol_initialize ( kjit, kjpindex, index, rest_id, & 438 SUBROUTINE hydrol_initialize ( ks, nvan, avan, mcr, & 439 mcs, mcfc, mcw, kjit, & 440 kjpindex, index, rest_id, & 453 441 njsc, soiltile, veget, veget_max, & 454 442 humrel, vegstress, drysoil_frac, & … … 461 449 !! 0. Variable and parameter declaration 462 450 !! 0.1 Input variables 451 452 453 !salma: added soil params as input variables 454 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ks !! Hydraulic conductivity at saturation (mm {-1}) 455 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: nvan !! Van Genuchten coeficients n (unitless) 456 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: avan !! Van Genuchten coeficients a (mm-1}) 457 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 458 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 459 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3}) 460 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3}) 461 !end salma: added soil params as input variables 462 463 463 INTEGER(i_std), INTENT(in) :: kjit !! Time step number 464 464 INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size … … 495 495 !! 0.4 Local variables 496 496 INTEGER(i_std) :: jsl 497 498 !salma: added soil params which depend on soil texture 499 REAL(r_std), DIMENSION (nscm) :: nvan_text 500 REAL(r_std), DIMENSION (nscm) :: avan_text 501 REAL(r_std), DIMENSION (nscm) :: mcr_text 502 REAL(r_std), DIMENSION (nscm) :: mcs_text 503 REAL(r_std), DIMENSION (nscm) :: ks_text 504 REAL(r_std), DIMENSION (nscm) :: pcent_text 505 REAL(r_std), DIMENSION (nscm) :: mcfc_text 506 REAL(r_std), DIMENSION (nscm) :: mcw_text 507 REAL(r_std), DIMENSION (nscm) :: mc_awet_text 508 REAL(r_std), DIMENSION (nscm) :: mc_adry_text 509 !end salma: added soil params which depend on soil texture 497 510 !_ ================================================================================================================================ 498 511 499 CALL hydrol_init (k jit, kjpindex, index, rest_id, veget_max, soiltile, &512 CALL hydrol_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc, kjit, kjpindex, index, rest_id, veget_max, soiltile, & 500 513 humrel, vegstress, snow, snow_age, snow_nobio, snow_nobio_age, qsintveg, & 501 514 snowdz, snowgrain, snowrho, snowtemp, snowheat, & 502 515 drysoil_frac, evap_bare_lim, evap_bare_lim_ns) 503 516 504 CALL hydrol_var_init (k jpindex, veget, veget_max, &517 CALL hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget, veget_max, & 505 518 soiltile, njsc, mx_eau_var, shumdiag_perma, & 506 519 drysoil_frac, qsintveg, mc_layh, mcl_layh) … … 557 570 !_ ================================================================================================================================ 558 571 559 SUBROUTINE hydrol_main (kjit, kjpindex, & 572 SUBROUTINE hydrol_main (ks, nvan, avan, mcr, mcs, mcfc, mcw, & 573 & kjit, kjpindex, & 560 574 & index, indexveg, indexsoil, indexlayer, indexnslm, & 561 575 & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, & … … 607 621 REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration 608 622 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: reinf_slope !! Slope coef 623 624 !salma: added input soil params: 625 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ks !! Hydraulic conductivity at saturation (mm {-1}) 626 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: nvan !! Van Genuchten coeficients n (unitless) 627 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: avan !! Van Genuchten coeficients a (mm-1}) 628 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 629 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 630 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3}) 631 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3}) 632 !end salma: added input soil params 633 609 634 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation 610 635 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_penm !! Soil Potential Evaporation Correction … … 672 697 673 698 !! 0.4 Local variables 674 699 !salma: added variable kfact_root_type 675 700 INTEGER(i_std) :: jst !! Index of soil tiles (unitless, 1-3) 676 701 INTEGER(i_std) :: jsl !! Index of soil layers (unitless) 677 702 INTEGER(i_std) :: ji, jv 703 CHARACTER(LEN=80) :: kfact_root_type !! read from run.def: when equal to 'cons', it indicates that ks does not increase 704 !! in the rootzone, ie, kfact_root=1; else, kfact_root defined as usual 678 705 REAL(r_std),DIMENSION (kjpindex) :: soilwet !! A temporary diagnostic of soil wetness 679 706 REAL(r_std),DIMENSION (kjpindex) :: snowdepth_diag !! Depth of snow layer containing default values, only for diagnostics … … 760 787 IF(soiltile(ji,jst) .GT. min_sechiba) THEN 761 788 kfact_root(ji,jsl,jst) = kfact_root(ji,jsl,jst) * & 762 & MAX((MAXVAL(ks_usda)/ks( njsc(ji)))**(- vegetmax_soil(ji,jv,jst)/2 * (humcste(jv)*zz(jsl)/mille - un)/deux), &789 & MAX((MAXVAL(ks_usda)/ks(ji))**(- vegetmax_soil(ji,jv,jst)/2 * (humcste(jv)*zz(jsl)/mille - un)/deux), & 763 790 un) 764 791 ENDIF … … 767 794 ENDDO 768 795 796 !salma: added config key of KFACT_ROOT_TYPE 797 !! KFACT_ROOT_TYPE = cons is used to impose that kfact_root = 1 in every soil layer, 798 !! so that ks does not increase in the rootzone; else, kfact_root defined as usual 799 !Config Key = KFACT_ROOT_TYPE 800 !Config Desc = keyword added for spmip exp1 and exp4 to get a constant ks over soil depth and rootzone 801 !Config If = spmip exp1 or exp4 802 !Config Def = var 803 !Config Help = can have two values: 'cons' or 'var'. If var then no changes, if cons then kfact_root=un 804 !Config Units = [mm/d] 805 CALL getin_p("KFACT_ROOT_TYPE",kfact_root_type) 806 807 IF (kfact_root_type=='cons') THEN 808 kfact_root(:,:,:) = un 809 ENDIF 810 !end salma: added config key 811 769 812 ! 770 813 !! 3.3 computes canopy ==>hydrol_canop … … 778 821 !! 3.5 computes soil hydrology ==>hydrol_soil 779 822 780 CALL hydrol_soil(k jpindex, veget_max, soiltile, njsc, reinf_slope, &823 CALL hydrol_soil(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, reinf_slope, & 781 824 transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 782 825 returnflow, reinfiltration, irrigation, & … … 869 912 ! mcs etc are identical in all layers (no normalization by vegtot to be comparable to mc) 870 913 DO jsl=1,nslm 871 land_mcs(:,jsl) = mcs(njsc(:)) 872 land_mcfc(:,jsl) = mcfc(njsc(:)) 873 land_mcw(:,jsl) = mcw(njsc(:)) 874 land_mcr(:,jsl) = mcr(njsc(:)) 914 !salma: replaced mcs(njsc(:)) by mcs(:) and same for other variables 915 land_mcs(:,jsl) = mcs(:) 916 land_mcfc(:,jsl) = mcfc(:) 917 land_mcw(:,jsl) = mcw(:) 918 land_mcr(:,jsl) = mcr(:) 875 919 ENDDO 876 920 CALL xios_orchidee_send_field("mcs",land_mcs) ! in m3/m3 … … 878 922 CALL xios_orchidee_send_field("mcw",land_mcw) ! in m3/m3 879 923 CALL xios_orchidee_send_field("mcr",land_mcr) ! in m3/m3 880 924 925 881 926 CALL xios_orchidee_send_field("water2infilt",water2infilt) 882 927 CALL xios_orchidee_send_field("mc",mc) … … 915 960 916 961 ! For the soil wetness above wilting point for CMIP6 (mrsow) 917 mrsow(:) = MAX( zero,humtot(:) - zmaxh*mille*mcw( njsc(:)) ) &918 / ( zmaxh*mille*( mcs( njsc(:)) - mcw(njsc(:)) ) )962 mrsow(:) = MAX( zero,humtot(:) - zmaxh*mille*mcw(:) ) & 963 / ( zmaxh*mille*( mcs(:) - mcw(:) ) ) 919 964 CALL xios_orchidee_send_field("mrsow",mrsow) 920 965 … … 1288 1333 !_ ================================================================================================================================ 1289 1334 !!_ hydrol_init 1290 1291 SUBROUTINE hydrol_init(kjit, kjpindex, index, rest_id, veget_max, soiltile, & 1335 !salma: added variables and njsc in arguments and input variables 1336 1337 SUBROUTINE hydrol_init(ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc,& 1338 kjit, kjpindex, index, rest_id, veget_max, soiltile, & 1292 1339 humrel, vegstress, snow, snow_age, snow_nobio, snow_nobio_age, qsintveg, & 1293 1340 snowdz, snowgrain, snowrho, snowtemp, snowheat, & … … 1298 1345 1299 1346 !! 0.1 Input variables 1300 1347 !salma introduced njsc variable 1348 INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 1301 1349 INTEGER(i_std), INTENT (in) :: kjit !! Time step number 1302 1350 INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size … … 1305 1353 REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Carte de vegetation max 1306 1354 REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless) 1307 1355 1308 1356 !! 0.2 Output variables 1357 1358 !salma: added input variables 1359 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ks !! Hydraulic conductivity at saturation (mm {-1}) 1360 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: nvan !! Van Genuchten coeficients n (unitless) 1361 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: avan !! Van Genuchten coeficients a (mm-1}) 1362 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 1363 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 1364 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3}) 1365 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3}) 1309 1366 1310 1367 REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Stress hydrique, relative humidity … … 1407 1464 !! 2.1 array allocation for soil texture 1408 1465 1409 ALLOCATE (nvan(nscm),stat=ier) 1410 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan','','') 1411 1412 ALLOCATE (avan(nscm),stat=ier) 1413 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan','','') 1414 1415 ALLOCATE (mcr(nscm),stat=ier) 1416 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcr','','') 1417 1418 ALLOCATE (mcs(nscm),stat=ier) 1419 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcs','','') 1420 1421 ALLOCATE (ks(nscm),stat=ier) 1422 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ks','','') 1423 1424 ALLOCATE (pcent(nscm),stat=ier) 1466 1467 ALLOCATE (pcent(kjpindex),stat=ier) 1425 1468 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable pcent','','') 1426 1427 ALLOCATE (mcfc(nscm),stat=ier) 1428 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcfc','','') 1429 1430 ALLOCATE (mcw(nscm),stat=ier) 1431 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcw','','') 1432 1433 ALLOCATE (mc_awet(nscm),stat=ier) 1469 1470 ALLOCATE (mc_awet(kjpindex),stat=ier) 1434 1471 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_awet','','') 1435 1472 1436 ALLOCATE (mc_adry( nscm),stat=ier)1473 ALLOCATE (mc_adry(kjpindex),stat=ier) 1437 1474 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_adry','','') 1438 1475 … … 1442 1479 CASE (3) 1443 1480 1444 nvan(:) = nvan_fao(:) 1445 avan(:) = avan_fao(:) 1446 mcr(:) = mcr_fao(:) 1447 mcs(:) = mcs_fao(:) 1448 ks(:) = ks_fao(:) 1449 pcent(:) = pcent_fao(:) 1450 mcfc(:) = mcf_fao(:) 1451 mcw(:) = mcw_fao(:) 1452 mc_awet(:) = mc_awet_fao(:) 1453 mc_adry(:) = mc_adry_fao(:) 1454 CASE (12) 1455 1456 nvan(:) = nvan_usda(:) 1457 avan(:) = avan_usda(:) 1458 mcr(:) = mcr_usda(:) 1459 mcs(:) = mcs_usda(:) 1460 ks(:) = ks_usda(:) 1461 pcent(:) = pcent_usda(:) 1462 mcfc(:) = mcf_usda(:) 1463 mcw(:) = mcw_usda(:) 1464 mc_awet(:) = mc_awet_usda(:) 1465 mc_adry(:) = mc_adry_usda(:) 1481 pcent(:) = pcent_fao(njsc(:)) 1482 mc_awet(:) = mc_awet_fao(njsc(:)) 1483 mc_adry(:) = mc_adry_fao(njsc(:)) 1484 CASE (13) 1485 1486 pcent(:) = pcent_usda(njsc(:)) 1487 mc_awet(:) = mc_awet_usda(njsc(:)) 1488 mc_adry(:) = mc_adry_usda(njsc(:)) 1466 1489 1467 1490 CASE DEFAULT … … 1474 1497 !! 2.3 Read in the run.def the parameters values defined by the user 1475 1498 1476 !Config Key = CWRR_N_VANGENUCHTEN 1477 !Config Desc = Van genuchten coefficient n 1478 !Config If = 1479 !Config Def = 1.89, 1.56, 1.31 1480 !Config Help = This parameter will be constant over the entire 1481 !Config simulated domain, thus independent from soil 1482 !Config texture. 1483 !Config Units = [-] 1484 CALL getin_p("CWRR_N_VANGENUCHTEN",nvan) 1485 1486 !! Check parameter value (correct range) 1487 IF ( ANY(nvan(:) <= zero) ) THEN 1488 CALL ipslerr_p(error_level, "hydrol_init.", & 1489 & "Wrong parameter value for CWRR_N_VANGENUCHTEN.", & 1490 & "This parameter should be positive. ", & 1491 & "Please, check parameter value in run.def. ") 1492 END IF 1493 1494 1495 !Config Key = CWRR_A_VANGENUCHTEN 1496 !Config Desc = Van genuchten coefficient a 1497 !Config If = 1498 !Config Def = 0.0075, 0.0036, 0.0019 1499 !Config Help = This parameter will be constant over the entire 1500 !Config simulated domain, thus independent from soil 1501 !Config texture. 1502 !Config Units = [1/mm] 1503 CALL getin_p("CWRR_A_VANGENUCHTEN",avan) 1504 1505 !! Check parameter value (correct range) 1506 IF ( ANY(avan(:) <= zero) ) THEN 1507 CALL ipslerr_p(error_level, "hydrol_init.", & 1508 & "Wrong parameter value for CWRR_A_VANGENUCHTEN.", & 1509 & "This parameter should be positive. ", & 1510 & "Please, check parameter value in run.def. ") 1511 END IF 1512 1513 1514 !Config Key = VWC_RESIDUAL 1515 !Config Desc = Residual soil water content 1516 !Config If = 1517 !Config Def = 0.065, 0.078, 0.095 1518 !Config Help = This parameter will be constant over the entire 1519 !Config simulated domain, thus independent from soil 1520 !Config texture. 1521 !Config Units = [m3/m3] 1522 CALL getin_p("VWC_RESIDUAL",mcr) 1523 1524 !! Check parameter value (correct range) 1525 IF ( ANY(mcr(:) < zero) .OR. ANY(mcr(:) > 1.) ) THEN 1526 CALL ipslerr_p(error_level, "hydrol_init.", & 1527 & "Wrong parameter value for VWC_RESIDUAL.", & 1528 & "This parameter is ranged between 0 and 1. ", & 1529 & "Please, check parameter value in run.def. ") 1530 END IF 1531 1532 1533 !Config Key = VWC_SAT 1534 !Config Desc = Saturated soil water content 1535 !Config If = 1536 !Config Def = 0.41, 0.43, 0.41 1537 !Config Help = This parameter will be constant over the entire 1538 !Config simulated domain, thus independent from soil 1539 !Config texture. 1540 !Config Units = [m3/m3] 1541 CALL getin_p("VWC_SAT",mcs) 1542 1543 !! Check parameter value (correct range) 1544 IF ( ANY(mcs(:) < zero) .OR. ANY(mcs(:) > 1.) .OR. ANY(mcs(:) <= mcr(:)) ) THEN 1545 CALL ipslerr_p(error_level, "hydrol_init.", & 1546 & "Wrong parameter value for VWC_SAT.", & 1547 & "This parameter should be greater than VWC_RESIDUAL and less than 1. ", & 1548 & "Please, check parameter value in run.def. ") 1549 END IF 1550 1551 1552 !Config Key = CWRR_KS 1553 !Config Desc = Hydraulic conductivity Saturation 1554 !Config If = 1555 !Config Def = 1060.8, 249.6, 62.4 1556 !Config Help = This parameter will be constant over the entire 1557 !Config simulated domain, thus independent from soil 1558 !Config texture. 1559 !Config Units = [mm/d] 1560 CALL getin_p("CWRR_KS",ks) 1561 1562 !! Check parameter value (correct range) 1563 IF ( ANY(ks(:) <= zero) ) THEN 1564 CALL ipslerr_p(error_level, "hydrol_init.", & 1565 & "Wrong parameter value for CWRR_KS.", & 1566 & "This parameter should be positive. ", & 1567 & "Please, check parameter value in run.def. ") 1568 END IF 1499 1500 1501 1502 1503 1504 1569 1505 1570 1506 … … 1587 1523 1588 1524 1589 !Config Key = VWC_FC 1590 !Config Desc = Volumetric water content field capacity 1591 !Config If = 1592 !Config Def = 0.32, 0.32, 0.32 1593 !Config Help = This parameter is independent from soil texture for 1594 !Config the time being. 1595 !Config Units = [m3/m3] 1596 CALL getin_p("VWC_FC",mcfc) 1597 1598 !! Check parameter value (correct range) 1599 IF ( ANY(mcfc(:) > mcs(:)) ) THEN 1600 CALL ipslerr_p(error_level, "hydrol_init.", & 1601 & "Wrong parameter value for VWC_FC.", & 1602 & "This parameter should be less than VWC_SAT. ", & 1603 & "Please, check parameter value in run.def. ") 1604 END IF 1605 1606 1607 !Config Key = VWC_WP 1608 !Config Desc = Volumetric water content Wilting pt 1609 !Config If = 1610 !Config Def = 0.10, 0.10, 0.10 1611 !Config Help = This parameter is independent from soil texture for 1612 !Config the time being. 1613 !Config Units = [m3/m3] 1614 CALL getin_p("VWC_WP",mcw) 1615 1616 !! Check parameter value (correct range) 1617 IF ( ANY(mcw(:) > mcfc(:)) .OR. ANY(mcw(:) < mcr(:)) ) THEN 1618 CALL ipslerr_p(error_level, "hydrol_init.", & 1619 & "Wrong parameter value for VWC_WP.", & 1620 & "This parameter should be greater or equal than VWC_RESIDUAL and less or equal than VWC_SAT.", & 1621 & "Please, check parameter value in run.def. ") 1622 END IF 1623 1525 1624 1526 1625 1527 !Config Key = VWC_MIN_FOR_WET_ALB … … 1745 1647 kk(:,:,:) = 276.48 1746 1648 1747 ALLOCATE (avan_mod_tab(nslm, nscm),stat=ier)1649 ALLOCATE (avan_mod_tab(nslm,kjpindex),stat=ier) 1748 1650 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan_mod_tab','','') 1749 1651 1750 ALLOCATE (nvan_mod_tab(nslm, nscm),stat=ier)1652 ALLOCATE (nvan_mod_tab(nslm,kjpindex),stat=ier) 1751 1653 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan_mod_tab','','') 1752 1654 … … 1932 1834 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact_root','','') 1933 1835 1934 ALLOCATE (kfact(nslm, nscm),stat=ier)1836 ALLOCATE (kfact(nslm, kjpindex),stat=ier) 1935 1837 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact','','') 1936 1838 … … 1944 1846 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dh','','') 1945 1847 1946 ALLOCATE (mc_lin(imin:imax, nscm),stat=ier)1848 ALLOCATE (mc_lin(imin:imax, kjpindex),stat=ier) 1947 1849 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_lin','','') 1948 1850 1949 ALLOCATE (k_lin(imin:imax, nslm, nscm),stat=ier)1851 ALLOCATE (k_lin(imin:imax, nslm, kjpindex),stat=ier) 1950 1852 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable k_lin','','') 1951 1853 1952 ALLOCATE (d_lin(imin:imax, nslm, nscm),stat=ier)1854 ALLOCATE (d_lin(imin:imax, nslm, kjpindex),stat=ier) 1953 1855 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable d_lin','','') 1954 1856 1955 ALLOCATE (a_lin(imin:imax, nslm, nscm),stat=ier)1857 ALLOCATE (a_lin(imin:imax, nslm, kjpindex),stat=ier) 1956 1858 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable a_lin','','') 1957 1859 1958 ALLOCATE (b_lin(imin:imax, nslm, nscm),stat=ier)1860 ALLOCATE (b_lin(imin:imax, nslm, kjpindex),stat=ier) 1959 1861 IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable b_lin','','') 1960 1862 … … 2451 2353 2452 2354 ! Allocation for soiltile related parameters 2453 IF ( ALLOCATED (nvan)) DEALLOCATE (nvan) 2454 IF ( ALLOCATED (avan)) DEALLOCATE (avan) 2455 IF ( ALLOCATED (mcr)) DEALLOCATE (mcr) 2456 IF ( ALLOCATED (mcs)) DEALLOCATE (mcs) 2457 IF ( ALLOCATED (ks)) DEALLOCATE (ks) 2355 2458 2356 IF ( ALLOCATED (pcent)) DEALLOCATE (pcent) 2459 IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc)2460 IF ( ALLOCATED (mcw)) DEALLOCATE (mcw)2461 2357 IF ( ALLOCATED (mc_awet)) DEALLOCATE (mc_awet) 2462 2358 IF ( ALLOCATED (mc_adry)) DEALLOCATE (mc_adry) … … 2811 2707 !_ hydrol_var_init 2812 2708 2813 SUBROUTINE hydrol_var_init (kjpindex, veget, veget_max, soiltile, njsc, & 2709 SUBROUTINE hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, & 2710 kjpindex, veget, veget_max, soiltile, njsc, & 2814 2711 mx_eau_var, shumdiag_perma, & 2815 2712 drysoil_frac, qsintveg, mc_layh, mcl_layh) … … 2820 2717 2821 2718 !! 0.1 Input variables 2719 !salma: added the following soil parameters 2720 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ks !! Hydraulic conductivity at saturation (mm {-1}) 2721 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: nvan !! Van Genuchten coeficients n (unitless) 2722 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: avan !! Van Genuchten coeficients a (mm-1}) 2723 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 2724 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 2725 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3}) 2726 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3}) 2822 2727 2823 2728 ! input scalar … … 2855 2760 REAL(r_std) :: nvan_mod !! VG parameter n modified from exponantial profile 2856 2761 !! (unitless) 2857 REAL(r_std), DIMENSION(nslm, nscm):: afact, nfact !! Multiplicative factor for decay of a and n with depth2762 REAL(r_std), DIMENSION(nslm,kjpindex) :: afact, nfact !! Multiplicative factor for decay of a and n with depth 2858 2763 !! (unitless) 2859 2764 ! parameters for "soil densification" with depth … … 2883 2788 REAL(r_std), DIMENSION (kjpindex,nslm) :: nvg !! VG param n modified with depth at each node 2884 2789 !! (unitless) 2885 INTEGER(i_std) :: jiref !! To identify the mc_lins where k_lin and d_lin 2790 !! need special treatment 2791 !salma: added local variable ii and jiref replaced with iiref 2792 INTEGER(i_std) :: ii 2793 INTEGER(i_std) :: iiref !! To identify the mc_lins where k_lin and d_lin 2886 2794 !! need special treatment 2887 2795 … … 3050 2958 END IF 3051 2959 2960 2961 3052 2962 !- 3053 2963 !! 3 Compute the profile for a and n 3054 2964 !- 3055 3056 ! For every soil texture 3057 DO jsc = 1, nscm 2965 !salma: for every grid cell 2966 DO ji = 1, kjpindex 3058 2967 DO jsl=1,nslm 3059 2968 ! PhD thesis of d'Orgeval, 2006, p81, Eq. 4.38; d'Orgeval et al. 2008, Eq. 2 3060 2969 ! Calibrated against Hapex-Sahel measurements 3061 kfact(jsl,jsc) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un) 3062 ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14 3063 3064 nfact(jsl,jsc) = ( kfact(jsl,jsc) )**nk_rel 3065 afact(jsl,jsc) = ( kfact(jsl,jsc) )**ak_rel 3066 ENDDO 3067 ENDDO 3068 3069 ! For every soil texture 3070 DO jsc = 1, nscm 2970 kfact(jsl,ji) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un) 2971 ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14 2972 2973 nfact(jsl,ji) = ( kfact(jsl,ji) )**nk_rel 2974 afact(jsl,ji) = ( kfact(jsl,ji) )**ak_rel 2975 ENDDO 2976 ENDDO 2977 !end salma 2978 2979 ! For every grid cell 2980 DO ji = 1, kjpindex 3071 2981 !- 3072 2982 !! 4 Compute the linearized values of k, a, b and d … … 3080 2990 3081 2991 ! We define 51 bounds for 50 bins of mc between mcr and mcs 3082 mc_lin(imin,j sc)=mcr(jsc)3083 mc_lin(imax,j sc)=mcs(jsc)3084 DO ji= imin+1, imax-1 ! ji=2,503085 mc_lin( ji,jsc) = mcr(jsc) + (ji-imin)*(mcs(jsc)-mcr(jsc))/(imax-imin)2992 mc_lin(imin,ji)=mcr(ji) 2993 mc_lin(imax,ji)=mcs(ji) 2994 DO ii= imin+1, imax-1 ! ii=2,50 2995 mc_lin(ii,ji) = mcr(ji) + (ii-imin)*(mcs(ji)-mcr(ji))/(imax-imin) 3086 2996 ENDDO 3087 2997 3088 2998 DO jsl = 1, nslm 3089 2999 ! From PhD thesis of d'Orgeval, 2006, p81, Eq. 4.42 3090 nvan_mod = n0 + (nvan(j sc)-n0) * nfact(jsl,jsc)3091 avan_mod = a0 + (avan(j sc)-a0) * afact(jsl,jsc)3000 nvan_mod = n0 + (nvan(ji)-n0) * nfact(jsl,ji) 3001 avan_mod = a0 + (avan(ji)-a0) * afact(jsl,ji) 3092 3002 m = un - un / nvan_mod 3093 3003 ! Creation of arrays for SP-MIP output by landpoint 3094 nvan_mod_tab(jsl,j sc) = nvan_mod3095 avan_mod_tab(jsl,j sc) = avan_mod3096 ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(j sc) * kfact(jsl,jsc)3097 DO ji = imax,imin,-13098 frac=MIN(un,(mc_lin( ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))3099 k_lin( ji,jsl,jsc) = ks(jsc) * kfact(jsl,jsc) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**23004 nvan_mod_tab(jsl,ji) = nvan_mod 3005 avan_mod_tab(jsl,ji) = avan_mod 3006 ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(ji) * kfact(jsl,ji) 3007 DO ii = imax,imin,-1 3008 frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji))) 3009 k_lin(ii,jsl,ji) = ks(ji) * kfact(jsl,ji) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2 3100 3010 ENDDO 3101 3011 3102 3012 ! k_lin should not be zero, nor too small 3103 ! We track jiref, the bin under which mc is too small and we may get zero k_lin 3104 ji=imax-1 3105 DO WHILE ((k_lin(ji,jsl,jsc) > 1.e-32) .and. (ji>0)) 3106 jiref=ji 3107 ji=ji-1 3013 ! We track iiref, the bin under which mc is too small and we may get zero k_lin 3014 !salma: ji replaced with ii and jiref replaced with iiref and jsc with ji 3015 ii=imax-1 3016 DO WHILE ((k_lin(ii,jsl,ji) > 1.e-32) .and. (ii>0)) 3017 iiref=ii 3018 ii=ii-1 3108 3019 ENDDO 3109 DO ji=jiref-1,imin,-13110 k_lin( ji,jsl,jsc)=k_lin(ji+1,jsl,jsc)/10.3020 DO ii=iiref-1,imin,-1 3021 k_lin(ii,jsl,ji)=k_lin(ii+1,jsl,ji)/10. 3111 3022 ENDDO 3112 3113 DO ji = imin,imax-1 ! ji=1,503023 3024 DO ii = imin,imax-1 ! ii=1,50 3114 3025 ! We deduce a_lin and b_lin based on continuity between segments k_lin = a_lin*mc-lin+b_lin 3115 a_lin( ji,jsl,jsc) = (k_lin(ji+1,jsl,jsc)-k_lin(ji,jsl,jsc)) / (mc_lin(ji+1,jsc)-mc_lin(ji,jsc))3116 b_lin( ji,jsl,jsc) = k_lin(ji,jsl,jsc) - a_lin(ji,jsl,jsc)*mc_lin(ji,jsc)3026 a_lin(ii,jsl,ji) = (k_lin(ii+1,jsl,ji)-k_lin(ii,jsl,ji)) / (mc_lin(ii+1,ji)-mc_lin(ii,ji)) 3027 b_lin(ii,jsl,ji) = k_lin(ii,jsl,ji) - a_lin(ii,jsl,ji)*mc_lin(ii,ji) 3117 3028 3118 3029 ! We calculate the d_lin for each mc bin, from Van Genuchten equation for D(theta) 3119 ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin 3120 IF ( ji.NE.imin .AND. ji.NE.imax-1) THEN3121 frac=MIN(un,(mc_lin( ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))3122 d_lin( ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) * &3123 ( (frac**(-un/m))/(mc_lin( ji,jsc)-mcr(jsc)) ) * &3030 ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin 3031 IF (ii.NE.imin .AND. ii.NE.imax-1) THEN 3032 frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji))) 3033 d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) * & 3034 ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) * & 3124 3035 ( frac**(-un/m) -un ) ** (-m) 3125 frac=MIN(un,(mc_lin( ji+1,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))3126 d_lin( ji+1,jsl,jsc) =(k_lin(ji+1,jsl,jsc) / (avan_mod*m*nvan_mod))*&3127 ( (frac**(-un/m))/(mc_lin( ji+1,jsc)-mcr(jsc)) ) * &3036 frac=MIN(un,(mc_lin(ii+1,ji)-mcr(ji))/(mcs(ji)-mcr(ji))) 3037 d_lin(ii+1,jsl,ji) =(k_lin(ii+1,jsl,ji) / (avan_mod*m*nvan_mod))*& 3038 ( (frac**(-un/m))/(mc_lin(ii+1,ji)-mcr(ji)) ) * & 3128 3039 ( frac**(-un/m) -un ) ** (-m) 3129 d_lin( ji,jsl,jsc) = undemi * (d_lin(ji,jsl,jsc)+d_lin(ji+1,jsl,jsc))3130 ELSE IF( ji.EQ.imax-1) THEN3131 d_lin( ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) * &3132 ( (frac**(-un/m))/(mc_lin( ji,jsc)-mcr(jsc)) ) * &3040 d_lin(ii,jsl,ji) = undemi * (d_lin(ii,jsl,ji)+d_lin(ii+1,jsl,ji)) 3041 ELSE IF(ii.EQ.imax-1) THEN 3042 d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) * & 3043 ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) * & 3133 3044 ( frac**(-un/m) -un ) ** (-m) 3134 3045 ENDIF 3135 ENDDO 3136 3137 ! Special case for ji=imin3138 d_lin(imin,jsl,j sc) = d_lin(imin+1,jsl,jsc)/1000.3046 ENDDO !Salma end loop over landpoints 3047 3048 ! Special case for ii=imin 3049 d_lin(imin,jsl,ji) = d_lin(imin+1,jsl,ji)/1000. 3139 3050 3140 3051 ! We adjust d_lin where k_lin was previously adjusted otherwise we might get non-monotonous variations 3141 3052 ! We don't want d_lin = zero 3142 DO ji=jiref-1,imin,-13143 d_lin( ji,jsl,jsc)=d_lin(ji+1,jsl,jsc)/10.3053 DO ii=iiref-1,imin,-1 3054 d_lin(ii,jsl,ji)=d_lin(ii+1,jsl,ji)/10. 3144 3055 ENDDO 3145 3056 3146 3057 ENDDO 3147 3058 ENDDO 3059 3148 3060 3149 3061 ! Output of alphavg and nvg at each node for SP-MIP 3150 3062 DO jsl = 1, nslm 3151 alphavg(:,jsl) = avan_mod_tab(jsl, njsc(:))*1000. ! from mm-1 to m-13152 nvg(:,jsl) = nvan_mod_tab(jsl, njsc(:))3063 alphavg(:,jsl) = avan_mod_tab(jsl,:)*1000. ! from mm-1 to m-1 3064 nvg(:,jsl) = nvan_mod_tab(jsl,:) 3153 3065 ENDDO 3154 3066 CALL xios_orchidee_send_field("alphavg",alphavg) ! in m-1 … … 3165 3077 3166 3078 mx_eau_var(:) = zero 3167 mx_eau_var(:) = zmaxh*mille*mcs( njsc(:))3168 3169 DO ji = 1,kjpindex 3079 mx_eau_var(:) = zmaxh*mille*mcs(:) 3080 3081 DO ji = 1,kjpindex 3170 3082 IF (vegtot(ji) .LE. zero) THEN 3171 3083 mx_eau_var(ji) = mx_eau_nobio*zmaxh … … 3181 3093 3182 3094 ! Loop on soiltiles to compute the variables (ji,jst) 3183 DO jst=1,nstm 3095 DO jst=1,nstm 3184 3096 DO ji = 1, kjpindex 3185 tmcs(ji,jst)=zmaxh* mille*mcs( njsc(ji))3186 tmcr(ji,jst)=zmaxh* mille*mcr( njsc(ji))3187 tmcfc(ji,jst)=zmaxh* mille*mcfc( njsc(ji))3188 tmcw(ji,jst)=zmaxh* mille*mcw( njsc(ji))3189 ENDDO 3190 ENDDO 3191 3097 tmcs(ji,jst)=zmaxh* mille*mcs(ji) 3098 tmcr(ji,jst)=zmaxh* mille*mcr(ji) 3099 tmcfc(ji,jst)=zmaxh* mille*mcfc(ji) 3100 tmcw(ji,jst)=zmaxh* mille*mcw(ji) 3101 ENDDO 3102 ENDDO 3103 3192 3104 ! The total soil moisture for each soiltile: 3193 3105 DO jst=1,nstm … … 3197 3109 ENDDO 3198 3110 3199 DO jst=1,nstm 3111 DO jst=1,nstm 3200 3112 DO jsl=2,nslm-1 3201 3113 DO ji=1,kjpindex … … 3206 3118 ENDDO 3207 3119 3208 DO jst=1,nstm 3120 DO jst=1,nstm 3209 3121 DO ji=1,kjpindex 3210 3122 tmc(ji,jst) = tmc(ji,jst) + dz(nslm) * (trois * mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit … … 3213 3125 END DO 3214 3126 3215 !JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty. 3127 !JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty. 3216 3128 ! ! If veget has been updated before restart (with LAND USE or DGVM), 3217 3129 ! ! tmc and mc must be modified with respect to humtot conservation. … … 3220 3132 ! The litter variables: 3221 3133 ! level 1 3222 DO jst=1,nstm 3134 DO jst=1,nstm 3223 3135 DO ji=1,kjpindex 3224 3136 tmc_litter(ji,jst) = dz(2) * (trois*mcl(ji,1,jst)+mcl(ji,2,jst))/huit 3225 tmc_litter_wilt(ji,jst) = dz(2) * mcw( njsc(ji)) / deux3226 tmc_litter_res(ji,jst) = dz(2) * mcr( njsc(ji)) / deux3227 tmc_litter_field(ji,jst) = dz(2) * mcfc( njsc(ji)) / deux3228 tmc_litter_sat(ji,jst) = dz(2) * mcs( njsc(ji)) / deux3229 tmc_litter_awet(ji,jst) = dz(2) * mc_awet( njsc(ji)) / deux3230 tmc_litter_adry(ji,jst) = dz(2) * mc_adry( njsc(ji)) / deux3137 tmc_litter_wilt(ji,jst) = dz(2) * mcw(ji) / deux 3138 tmc_litter_res(ji,jst) = dz(2) * mcr(ji) / deux 3139 tmc_litter_field(ji,jst) = dz(2) * mcfc(ji) / deux 3140 tmc_litter_sat(ji,jst) = dz(2) * mcs(ji) / deux 3141 tmc_litter_awet(ji,jst) = dz(2) * mc_awet(ji) / deux 3142 tmc_litter_adry(ji,jst) = dz(2) * mc_adry(ji) / deux 3231 3143 ENDDO 3232 3144 END DO 3233 3145 ! sum from level 2 to 4 3234 DO jst=1,nstm 3146 DO jst=1,nstm 3235 3147 DO jsl=2,4 3236 3148 DO ji=1,kjpindex 3237 tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * & 3149 tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * & 3238 3150 & ( trois*mcl(ji,jsl,jst) + mcl(ji,jsl-1,jst))/huit & 3239 3151 & + dz(jsl+1)*(trois*mcl(ji,jsl,jst) + mcl(ji,jsl+1,jst))/huit 3240 3152 tmc_litter_wilt(ji,jst) = tmc_litter_wilt(ji,jst) + & 3241 &(dz(jsl)+ dz(jsl+1))*& 3242 & mcw( njsc(ji))/deux3153 &(dz(jsl)+ dz(jsl+1))*& 3154 & mcw(ji)/deux 3243 3155 tmc_litter_res(ji,jst) = tmc_litter_res(ji,jst) + & 3244 &(dz(jsl)+ dz(jsl+1))*& 3245 & mcr( njsc(ji))/deux3156 &(dz(jsl)+ dz(jsl+1))*& 3157 & mcr(ji)/deux 3246 3158 tmc_litter_sat(ji,jst) = tmc_litter_sat(ji,jst) + & 3247 &(dz(jsl)+ dz(jsl+1))* & 3248 & mcs( njsc(ji))/deux3159 &(dz(jsl)+ dz(jsl+1))* & 3160 & mcs(ji)/deux 3249 3161 tmc_litter_field(ji,jst) = tmc_litter_field(ji,jst) + & 3250 & (dz(jsl)+ dz(jsl+1))* & 3251 & mcfc( njsc(ji))/deux3162 & (dz(jsl)+ dz(jsl+1))* & 3163 & mcfc(ji)/deux 3252 3164 tmc_litter_awet(ji,jst) = tmc_litter_awet(ji,jst) + & 3253 &(dz(jsl)+ dz(jsl+1))* & 3254 & mc_awet( njsc(ji))/deux3165 &(dz(jsl)+ dz(jsl+1))* & 3166 & mc_awet(ji)/deux 3255 3167 tmc_litter_adry(ji,jst) = tmc_litter_adry(ji,jst) + & 3256 & (dz(jsl)+ dz(jsl+1))* & 3257 & mc_adry( njsc(ji))/deux3168 & (dz(jsl)+ dz(jsl+1))* & 3169 & mc_adry(ji)/deux 3258 3170 END DO 3259 3171 END DO … … 3261 3173 3262 3174 3263 DO jst=1,nstm 3175 DO jst=1,nstm 3264 3176 DO ji=1,kjpindex 3265 3177 ! here we set that humrelv=0 in PFT1 3266 3178 humrelv(ji,1,jst) = zero 3267 3179 ENDDO 3268 3180 END DO … … 3270 3182 3271 3183 ! Calculate shumdiag_perma for thermosoil 3272 ! Use resdist instead of soiltile because we here need to have 3184 ! Use resdist instead of soiltile because we here need to have 3273 3185 ! shumdiag_perma at the value from previous time step. 3274 3186 ! Here, soilmoist is only used as a temporary variable to calculate shumdiag_perma … … 3290 3202 ENDDO 3291 3203 DO ji=1,kjpindex 3292 soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average 3293 ENDDO 3294 3204 soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average 3205 ENDDO 3206 3295 3207 ! -- shumdiag_perma for restart 3296 3208 ! For consistency with hydrol_soil, we want to calculate a grid-cell average 3297 3209 DO jsl = 1, nslm 3298 DO ji=1,kjpindex 3299 shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs( njsc(ji)))3300 shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero) 3301 ENDDO 3302 ENDDO 3303 3210 DO ji=1,kjpindex 3211 shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(ji)) 3212 shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero) 3213 ENDDO 3214 ENDDO 3215 3304 3216 ! Calculate drysoil_frac if it was not found in the restart file 3305 3217 ! For simplicity, we set drysoil_frac to 0.5 in this case … … 3310 3222 END IF 3311 3223 3312 !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in 3313 !! thermosoil for the thermal conductivity. 3224 !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in 3225 !! thermosoil for the thermal conductivity. 3314 3226 ! These values are only used in thermosoil_init in absence of a restart file 3227 3315 3228 mc_layh(:,:) = zero 3316 3229 mcl_layh(:,:) = zero 3230 3317 3231 DO jst=1,nstm 3318 3232 DO jsl=1,nslm … … 3668 3582 3669 3583 END SUBROUTINE hydrol_flood 3670 3671 3584 3672 3585 !! ================================================================================================================================ … … 3734 3647 !_ ================================================================================================================================ 3735 3648 !_ hydrol_soil 3736 3737 SUBROUTINE hydrol_soil (kjpindex, veget_max, soiltile, njsc, reinf_slope, &3649 SUBROUTINE hydrol_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, & 3650 kjpindex, veget_max, soiltile, njsc, reinf_slope, & 3738 3651 & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 3739 3652 & returnflow, reinfiltration, irrigation, & … … 3748 3661 3749 3662 !! 0.1 Input variables 3663 !salma added soil params to arguments and input variables 3664 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ks !! Hydraulic conductivity at saturation (mm {-1}) 3665 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: nvan !! Van Genuchten coeficients n (unitless) 3666 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: avan !! Van Genuchten coeficients a (mm-1}) 3667 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 3668 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 3669 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3}) 3670 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3}) 3671 3750 3672 3751 3673 INTEGER(i_std), INTENT(in) :: kjpindex … … 3807 3729 !! averaged over the mesh (for thermosoil) 3808 3730 !! @tex $(m^{3} m^{-3})$ @endtex 3809 3810 3731 !! 0.3 Modified variables 3811 3732 … … 3953 3874 ! These values will be kept till the end of the prognostic loop 3954 3875 DO jst=1,nstm 3955 CALL hydrol_soil_froz( kjpindex,jst,njsc)3876 CALL hydrol_soil_froz(nvan, avan, mcr, mcs, kjpindex,jst,njsc) 3956 3877 ENDDO 3957 3878 … … 4080 4001 !! K and D are computed for the profile of mc before infiltration 4081 4002 !! They depend on the fraction of soil ice, given by profil_froz_hydro_ns 4082 CALL hydrol_soil_coef( kjpindex,jst,njsc)4003 CALL hydrol_soil_coef(mcr, mcs, kjpindex,jst,njsc) 4083 4004 4084 4005 !! Infiltration and surface runoff are computed … … 4086 4007 !! The conductivity comes from hydrol_soil_coef and relates to the liquid phase only 4087 4008 ! This seems consistent with ok_freeze 4088 CALL hydrol_soil_infilt(k jpindex, jst, njsc, flux_infilt, qinfilt_ns, ru_infilt_ns, &4009 CALL hydrol_soil_infilt(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, jst, njsc, flux_infilt, qinfilt_ns, ru_infilt_ns, & 4089 4010 check_infilt_ns) 4090 4011 ru_ns(:,jst) = ru_infilt_ns(:,jst) … … 4112 4033 !! 2.1 K and D are recomputed after infiltration 4113 4034 !! They depend on the fraction of soil ice, still given by profil_froz_hydro_ns 4114 CALL hydrol_soil_coef( kjpindex,jst,njsc)4035 CALL hydrol_soil_coef(mcr, mcs,kjpindex,jst,njsc) 4115 4036 4116 4037 !! 2.2 Set the tridiagonal matrix coefficients for the diffusion/redistribution scheme … … 4121 4042 DO jsl = 1, nslm 4122 4043 DO ji =1, kjpindex 4123 mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr( njsc(ji)) + &4124 (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr( njsc(ji))) )4044 mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + & 4045 (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) ) 4125 4046 ! we always have mcl<=mc 4126 4047 ! if mc>mcr, then mcl>mcr; if mc=mcr,mcl=mcr; if mc<mcr, then mcl<mcr … … 4143 4064 DO ji = 1, kjpindex 4144 4065 IF (soiltile(ji,jst) .GT. zero) THEN 4145 mvg = un - un / nvan_mod_tab(jsl, njsc(ji))4146 avg = avan_mod_tab(jsl, njsc(ji))*1000. ! to convert in m-14147 mc_ratio = MAX( 10.**(-14*mvg), (mcl(ji,jsl,jst) - mcr( njsc(ji)))/(mcs(njsc(ji)) - mcr(njsc(ji))) )**(-un/mvg)4148 psi = - MAX(zero,(mc_ratio - un))**(un/nvan_mod_tab(jsl, njsc(ji))) / avg ! in m4066 mvg = un - un / nvan_mod_tab(jsl,ji) 4067 avg = avan_mod_tab(jsl,ji)*1000. ! to convert in m-1 4068 mc_ratio = MAX( 10.**(-14*mvg), (mcl(ji,jsl,jst) - mcr(ji))/(mcs(ji) - mcr(ji)) )**(-un/mvg) 4069 psi = - MAX(zero,(mc_ratio - un))**(un/nvan_mod_tab(jsl,ji)) / avg ! in m 4149 4070 psi_moy(ji,jsl) = psi_moy(ji,jsl) + soiltile(ji,jst) * psi ! average across soil tiles 4150 4071 ENDIF … … 4274 4195 DO ji =1, kjpindex 4275 4196 mc(ji,jsl,jst) = MAX( mcl(ji,jsl,jst), mcl(ji,jsl,jst) + & 4276 profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr( njsc(ji))) )4197 profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(ji)) ) 4277 4198 ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we get mc=mcl 4278 4199 ENDDO … … 4286 4207 dr_corr_ns(:,jst) = zero 4287 4208 ru_corr_ns(:,jst) = zero 4288 call hydrol_soil_smooth_over_mcs2( kjpindex, jst, njsc, is_over_mcs, ru_corr_ns, check_over_ns)4209 call hydrol_soil_smooth_over_mcs2(mcs, kjpindex, jst, njsc, is_over_mcs, ru_corr_ns, check_over_ns) 4289 4210 4290 4211 ! In absence of freezing, if F is large enough, the correction of oversaturation is sent to drainage … … 4330 4251 dmc(ji,jsl) = zero 4331 4252 IF ( ( zz(jsl) >= zwt_force(ji,jst)*mille ) ) THEN 4332 dmc(ji,jsl) = mcs( njsc(ji)) - mc(ji,jsl,jst) ! addition to reach mcs (m3/m3) = positive value4333 mc(ji,jsl,jst) = mcs( njsc(ji))4253 dmc(ji,jsl) = mcs(ji) - mc(ji,jsl,jst) ! addition to reach mcs (m3/m3) = positive value 4254 mc(ji,jsl,jst) = mcs(ji) 4334 4255 ENDIF 4335 4256 ENDDO … … 4366 4287 wtd_ns(ji,jst) = undef_sechiba ! in meters 4367 4288 jsl=nslm 4368 DO WHILE ( (mc(ji,jsl,jst) .EQ. mcs( njsc(ji))) .AND. (jsl > 1) )4289 DO WHILE ( (mc(ji,jsl,jst) .EQ. mcs(ji)) .AND. (jsl > 1) ) 4369 4290 wtd_ns(ji,jst) = zz(jsl)/mille ! in meters 4370 4291 jsl=jsl-1 … … 4375 4296 ! This routine does not change tmc but decides where we should turn off ET to prevent further mc decrease 4376 4297 ! Like above, the tests are made on total mc, compared to mcr 4377 CALL hydrol_soil_smooth_under_mcr( kjpindex, jst, njsc, is_under_mcr, check_under_ns)4298 CALL hydrol_soil_smooth_under_mcr(mcr, kjpindex, jst, njsc, is_under_mcr, check_under_ns) 4378 4299 4379 4300 !! 4. At the end of the prognostic calculations, we recompute important moisture variables … … 4399 4320 DO jsl = 1, nslm 4400 4321 DO ji =1, kjpindex 4401 mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr( njsc(ji)) + &4402 (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr( njsc(ji))) )4322 mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + & 4323 (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) ) 4403 4324 ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we keep mcl=mc 4404 4325 ENDDO … … 4451 4372 sm(:,1) = dz(2) * (trois*mcl(:,1,jst) + mcl(:,2,jst))/huit 4452 4373 smt(:,1) = dz(2) * (trois*mc(:,1,jst) + mc(:,2,jst))/huit 4453 smw(:,1) = dz(2) * (quatre*mcw( njsc(:)))/huit4454 smf(:,1) = dz(2) * (quatre*mcfc( njsc(:)))/huit4455 sms(:,1) = dz(2) * (quatre*mcs( njsc(:)))/huit4374 smw(:,1) = dz(2) * (quatre*mcw(:))/huit 4375 smf(:,1) = dz(2) * (quatre*mcfc(:))/huit 4376 sms(:,1) = dz(2) * (quatre*mcs(:))/huit 4456 4377 DO jsl = 2,nslm-1 4457 4378 sm(:,jsl) = dz(jsl) * (trois*mcl(:,jsl,jst)+mcl(:,jsl-1,jst))/huit & … … 4459 4380 smt(:,jsl) = dz(jsl) * (trois*mc(:,jsl,jst)+mc(:,jsl-1,jst))/huit & 4460 4381 + dz(jsl+1) * (trois*mc(:,jsl,jst)+mc(:,jsl+1,jst))/huit 4461 smw(:,jsl) = dz(jsl) * ( quatre*mcw( njsc(:)) )/huit &4462 + dz(jsl+1) * ( quatre*mcw( njsc(:)) )/huit4463 smf(:,jsl) = dz(jsl) * ( quatre*mcfc( njsc(:)) )/huit &4464 + dz(jsl+1) * ( quatre*mcfc( njsc(:)) )/huit4465 sms(:,jsl) = dz(jsl) * ( quatre*mcs( njsc(:)) )/huit &4466 + dz(jsl+1) * ( quatre*mcs( njsc(:)) )/huit4382 smw(:,jsl) = dz(jsl) * ( quatre*mcw(:) )/huit & 4383 + dz(jsl+1) * ( quatre*mcw(:) )/huit 4384 smf(:,jsl) = dz(jsl) * ( quatre*mcfc(:) )/huit & 4385 + dz(jsl+1) * ( quatre*mcfc(:) )/huit 4386 sms(:,jsl) = dz(jsl) * ( quatre*mcs(:) )/huit & 4387 + dz(jsl+1) * ( quatre*mcs(:) )/huit 4467 4388 ENDDO 4468 4389 sm(:,nslm) = dz(nslm) * (trois*mcl(:,nslm,jst) + mcl(:,nslm-1,jst))/huit 4469 4390 smt(:,nslm) = dz(nslm) * (trois*mc(:,nslm,jst) + mc(:,nslm-1,jst))/huit 4470 smw(:,nslm) = dz(nslm) * (quatre*mcw( njsc(:)))/huit4471 smf(:,nslm) = dz(nslm) * (quatre*mcfc( njsc(:)))/huit4472 sms(:,nslm) = dz(nslm) * (quatre*mcs( njsc(:)))/huit4391 smw(:,nslm) = dz(nslm) * (quatre*mcw(:))/huit 4392 smf(:,nslm) = dz(nslm) * (quatre*mcfc(:))/huit 4393 sms(:,nslm) = dz(nslm) * (quatre*mcs(:))/huit 4473 4394 ! sm_nostress = soil moisture of each layer at which us reaches 1, here at the middle of [smw,smf] 4474 4395 DO jsl = 1,nslm 4475 sm_nostress(:,jsl) = smw(:,jsl) + pcent( njsc(:)) * (smf(:,jsl)-smw(:,jsl))4396 sm_nostress(:,jsl) = smw(:,jsl) + pcent(:) * (smf(:,jsl)-smw(:,jsl)) 4476 4397 END DO 4477 4398 … … 4642 4563 !! The multiplication by vegtot creates grid-cell average values 4643 4564 ! *** To be checked for consistency with the use of nobio properties in thermosoil 4565 4644 4566 DO jsl=1,nslm 4645 4567 DO ji=1,kjpindex … … 4661 4583 DO jsl = 1, nslm 4662 4584 ksat(:,jsl) = ksat(:,jsl) + soiltile(:,jst) * & 4663 ( ks( njsc(:)) * kfact(jsl,njsc(:)) * kfact_root(:,jsl,jst) )4585 ( ks(:) * kfact(jsl,:) * kfact_root(:,jsl,jst) ) 4664 4586 ENDDO 4665 4587 … … 4671 4593 4672 4594 !! 7. Summing 3d variables into 2d variables 4673 CALL hydrol_diag_soil (k jpindex, veget_max, soiltile, njsc, runoff, drainage, &4595 CALL hydrol_diag_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, runoff, drainage, & 4674 4596 & evapot, vevapnu, returnflow, reinfiltration, irrigation, & 4675 4597 & shumdiag,shumdiag_perma, k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,tot_melt) … … 4762 4684 !! 8.2 Since we estimate bare soile evap for the next time step, we update profil_froz_hydro and mcl 4763 4685 ! (effect of mc only, the change in temp_hydro is neglected) 4764 IF ( ok_freeze_cwrr ) CALL hydrol_soil_froz( kjpindex,jst,njsc)4686 IF ( ok_freeze_cwrr ) CALL hydrol_soil_froz(nvan, avan, mcr, mcs,kjpindex,jst,njsc) 4765 4687 DO jsl = 1, nslm 4766 4688 DO ji =1, kjpindex 4767 mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr( njsc(ji)) + &4768 (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr( njsc(ji))) )4689 mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + & 4690 (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) ) 4769 4691 ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we keep mcl=mc 4770 4692 ENDDO … … 4772 4694 4773 4695 !! 8.3 K and D are recomputed for the updated profile of mc/mcl 4774 CALL hydrol_soil_coef( kjpindex,jst,njsc)4696 CALL hydrol_soil_coef(mcr, mcs, kjpindex,jst,njsc) 4775 4697 4776 4698 !! 8.4 Set the tridiagonal matrix coefficients for the diffusion/redistribution scheme … … 4811 4733 4812 4734 flux_top(ji) = zero 4735 r_soil_ns(ji,jst) = zero 4813 4736 4814 4737 ENDIF … … 4863 4786 DO ji = 1, kjpindex 4864 4787 ! by construction, mc and mcl are always on the same side of mcr, so we can use mcl here 4865 resolv(ji) = (mcl(ji,1,jst).LT.(mcr( njsc(ji))).AND.flux_top(ji).GT.min_sechiba)4788 resolv(ji) = (mcl(ji,1,jst).LT.(mcr(ji)).AND.flux_top(ji).GT.min_sechiba) 4866 4789 ENDDO 4867 4790 !! Reset the coefficient for diffusion (tridiag is only solved if resolv(ji) = .TRUE.)O … … 4879 4802 tmat(ji,1,2) = un 4880 4803 tmat(ji,1,3) = zero 4881 rhs(ji,1) = mcr( njsc(ji))4804 rhs(ji,1) = mcr(ji) 4882 4805 ENDDO 4883 4806 … … 4898 4821 DO ji =1, kjpindex 4899 4822 mc(ji,jsl,jst) = MAX( mcl(ji,jsl,jst), mcl(ji,jsl,jst) + & 4900 profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr( njsc(ji))) )4823 profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(ji)) ) 4901 4824 ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we get mc=mcl 4902 4825 ENDDO … … 5050 4973 !_ hydrol_soil_infilt 5051 4974 5052 SUBROUTINE hydrol_soil_infilt(k jpindex, ins, njsc, flux_infilt, qinfilt_ns, ru_infilt, check)4975 SUBROUTINE hydrol_soil_infilt(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, ins, njsc, flux_infilt, qinfilt_ns, ru_infilt, check) 5053 4976 5054 4977 !! 0. Variable and parameter declaration … … 5057 4980 5058 4981 ! GLOBAL (in or inout) 4982 ! salma added input soil hydraulic parameters 4983 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ks !! Hydraulic conductivity at saturation (mm {-1}) 4984 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: nvan !! Van Genuchten coeficients n (unitless) 4985 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: avan !! Van Genuchten coeficients a (mm-1}) 4986 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 4987 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 4988 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3}) 4989 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3}) 4990 4991 5059 4992 INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size 5060 4993 INTEGER(i_std), INTENT(in) :: ins … … 5108 5041 DO ji = 1, kjpindex 5109 5042 !! First we fill up the first layer (about 1mm) without any resistance and quasi-immediately 5110 wat_inf_pot(ji) = MAX((mcs( njsc(ji))-mc(ji,1,ins)) * dz(2) / deux, zero)5043 wat_inf_pot(ji) = MAX((mcs(ji)-mc(ji,1,ins)) * dz(2) / deux, zero) 5111 5044 wat_inf(ji) = MIN(wat_inf_pot(ji), flux_infilt(ji)) 5112 5045 mc(ji,1,ins) = mc(ji,1,ins) + wat_inf(ji) * deux / dz(2) … … 5127 5060 ! This is computed by an simple arithmetic average because 5128 5061 ! the time step (30min) is not appropriate for a geometric average (advised by Haverkamp and Vauclin) 5129 k_m = (k(ji,jsl) + ks( njsc(ji))*kfact(jsl-1,njsc(ji))*kfact_root(ji,jsl,ins)) / deux5062 k_m = (k(ji,jsl) + ks(ji)*kfact(jsl-1,ji)*kfact_root(ji,jsl,ins)) / deux 5130 5063 5131 5064 IF (ok_freeze_cwrr) THEN … … 5141 5074 5142 5075 !! From which we deduce the time it takes to fill up the layer or to end the time step... 5143 wat_inf_pot(ji) = MAX((mcs( njsc(ji))-mc(ji,jsl,ins)) * (dz(jsl) + dz(jsl+1)) / deux, zero)5076 wat_inf_pot(ji) = MAX((mcs(ji)-mc(ji,jsl,ins)) * (dz(jsl) + dz(jsl+1)) / deux, zero) 5144 5077 IF ( infilt_tmp(ji) > min_sechiba) THEN 5145 5078 dt_inf(ji) = MIN(wat_inf_pot(ji)/infilt_tmp(ji), dt_tmp(ji)) … … 5223 5156 !_ hydrol_soil_smooth_under_mcr 5224 5157 5225 SUBROUTINE hydrol_soil_smooth_under_mcr( kjpindex, ins, njsc, is_under_mcr, check)5158 SUBROUTINE hydrol_soil_smooth_under_mcr(mcr, kjpindex, ins, njsc, is_under_mcr, check) 5226 5159 5227 5160 !- arguments … … 5231 5164 !! 0.1 Input variables 5232 5165 5166 !salma: added parameter mcr 5167 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 5233 5168 INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size 5234 5169 INTEGER(i_std), INTENT(in) :: ins !! Soiltile index (1-nstm, unitless) … … 5268 5203 DO jsl = 1,nslm-2 5269 5204 DO ji=1, kjpindex 5270 excess = MAX(mcr( njsc(ji))-mc(ji,jsl,ins),zero)5205 excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero) 5271 5206 mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess 5272 5207 mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) - excess * & … … 5277 5212 jsl = nslm-1 5278 5213 DO ji=1, kjpindex 5279 excess = MAX(mcr( njsc(ji))-mc(ji,jsl,ins),zero)5214 excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero) 5280 5215 mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess 5281 5216 mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) - excess * & … … 5285 5220 jsl = nslm 5286 5221 DO ji=1, kjpindex 5287 excess = MAX(mcr( njsc(ji))-mc(ji,jsl,ins),zero)5222 excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero) 5288 5223 mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess 5289 5224 mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) - excess * & … … 5294 5229 DO jsl = nslm-1,2,-1 5295 5230 DO ji=1, kjpindex 5296 excess = MAX(mcr( njsc(ji))-mc(ji,jsl,ins),zero)5231 excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero) 5297 5232 mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess 5298 5233 mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) - excess * & … … 5304 5239 ! excess > 0 5305 5240 DO ji=1, kjpindex 5306 excessji(ji) = mask_soiltile(ji,ins) * MAX(mcr( njsc(ji))-mc(ji,1,ins),zero)5241 excessji(ji) = mask_soiltile(ji,ins) * MAX(mcr(ji)-mc(ji,1,ins),zero) 5307 5242 ENDDO 5308 5243 DO ji=1, kjpindex … … 5380 5315 !_ hydrol_soil_smooth_over_mcs 5381 5316 5382 SUBROUTINE hydrol_soil_smooth_over_mcs( kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)5317 SUBROUTINE hydrol_soil_smooth_over_mcs(mcs ,kjpindex, ins, njsc, is_over_mcs, rudr_corr, check) 5383 5318 5384 5319 !- arguments … … 5386 5321 !! 0. Variable and parameter declaration 5387 5322 5388 !! 0.1 Input variables 5389 5323 !! 0.1 Input variables 5324 !salma: added mcs 5325 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 5390 5326 INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size 5391 5327 INTEGER(i_std), INTENT(in) :: ins !! Soiltile index (1-nstm, unitless) … … 5425 5361 DO jsl = 1, nslm-2 5426 5362 DO ji=1, kjpindex 5427 excess = MAX(mc(ji,jsl,ins)-mcs( njsc(ji)),zero)5363 excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero) 5428 5364 mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess 5429 5365 mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) + excess * & … … 5434 5370 jsl = nslm-1 5435 5371 DO ji=1, kjpindex 5436 excess = MAX(mc(ji,jsl,ins)-mcs( njsc(ji)),zero)5372 excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero) 5437 5373 mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess 5438 5374 mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) + excess * & … … 5442 5378 jsl = nslm 5443 5379 DO ji=1, kjpindex 5444 excess = MAX(mc(ji,jsl,ins)-mcs( njsc(ji)),zero)5380 excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero) 5445 5381 mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess 5446 5382 mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) + excess * & … … 5451 5387 DO jsl = nslm-1,2,-1 5452 5388 DO ji=1, kjpindex 5453 excess = MAX(mc(ji,jsl,ins)-mcs( njsc(ji)),zero)5389 excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero) 5454 5390 mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess 5455 5391 mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) + excess * & … … 5461 5397 5462 5398 DO ji=1, kjpindex 5463 excess = mask_soiltile(ji,ins) * MAX(mc(ji,1,ins)-mcs( njsc(ji)),zero)5399 excess = mask_soiltile(ji,ins) * MAX(mc(ji,1,ins)-mcs(ji),zero) 5464 5400 mc(ji,1,ins) = mc(ji,1,ins) - excess ! then mc(1)=mcs 5465 5401 rudr_corr(ji,ins) = rudr_corr(ji,ins) + excess * dz(2) / deux … … 5512 5448 !_ hydrol_soil_smooth_over_mcs2 5513 5449 5514 SUBROUTINE hydrol_soil_smooth_over_mcs2( kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)5450 SUBROUTINE hydrol_soil_smooth_over_mcs2(mcs, kjpindex, ins, njsc, is_over_mcs, rudr_corr, check) 5515 5451 5516 5452 !- arguments … … 5519 5455 5520 5456 !! 0.1 Input variables 5521 5457 !salma: added mcs 5458 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 5522 5459 INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size 5523 5460 INTEGER(i_std), INTENT(in) :: ins !! Soiltile index (1-nstm, unitless) … … 5562 5499 DO jsl = 1, nslm 5563 5500 DO ji=1, kjpindex 5564 excess(ji,jsl) = MAX(mc(ji,jsl,ins)-mcs( njsc(ji)),zero) ! >=05501 excess(ji,jsl) = MAX(mc(ji,jsl,ins)-mcs(ji),zero) ! >=0 5565 5502 mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess(ji,jsl) ! here mc either does not change or decreases 5566 5503 ENDDO … … 5796 5733 !_ ================================================================================================================================ 5797 5734 !_ hydrol_soil_coef 5798 5799 SUBROUTINE hydrol_soil_coef( kjpindex,ins,njsc)5735 5736 SUBROUTINE hydrol_soil_coef(mcr, mcs, kjpindex,ins,njsc) 5800 5737 5801 5738 IMPLICIT NONE … … 5804 5741 5805 5742 !! 0.1 Input variables 5806 5743 !salma: added mcr and mcs 5744 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 5745 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 5807 5746 INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size 5808 5747 INTEGER(i_std), INTENT(in) :: ins !! Index of soil type … … 5836 5775 ! It corresponds to a liquid mc, but the expression is different from mcl in hydrol_soil, 5837 5776 ! to ensure that we get the a, b, d of the first bin when mcl<mcr 5838 mc_used = mcr( njsc(ji))+x*MAX((mc(ji,jsl, ins)-mcr(njsc(ji))),zero)5777 mc_used = mcr(ji)+x*MAX((mc(ji,jsl, ins)-mcr(ji)),zero) 5839 5778 ! 5840 5779 ! calcul de k based on mc_liq 5841 5780 ! 5842 i= MAX(imin, MIN(imax-1, INT(imin +(imax-imin)*(mc_used-mcr( njsc(ji)))/(mcs(njsc(ji))-mcr(njsc(ji))))))5843 a(ji,jsl) = a_lin(i,jsl, njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d5844 b(ji,jsl) = b_lin(i,jsl, njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d5845 d(ji,jsl) = d_lin(i,jsl, njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d5846 k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl, njsc(ji)), &5847 a_lin(i,jsl, njsc(ji)) * mc_used + b_lin(i,jsl,njsc(ji))) ! in mm/d5781 i= MAX(imin, MIN(imax-1, INT(imin +(imax-imin)*(mc_used-mcr(ji))/(mcs(ji)-mcr(ji))))) 5782 a(ji,jsl) = a_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d 5783 b(ji,jsl) = b_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d 5784 d(ji,jsl) = d_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm^2/d 5785 k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,ji), & 5786 a_lin(i,jsl,ji) * mc_used + b_lin(i,jsl,ji)) ! in mm/d 5848 5787 ENDDO ! loop on grid 5849 5788 ENDDO … … 5855 5794 5856 5795 ! it is impossible to consider a mc<mcr for the binning 5857 mc_ratio = MAX(mc(ji,jsl,ins)-mcr( njsc(ji)), zero)/(mcs(njsc(ji))-mcr(njsc(ji)))5796 mc_ratio = MAX(mc(ji,jsl,ins)-mcr(ji), zero)/(mcs(ji)-mcr(ji)) 5858 5797 5859 5798 i= MAX(MIN(INT((imax-imin)*mc_ratio)+imin , imax-1), imin) 5860 a(ji,jsl) = a_lin(i,jsl, njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d5861 b(ji,jsl) = b_lin(i,jsl, njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d5862 d(ji,jsl) = d_lin(i,jsl, njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d5863 k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl, njsc(ji)), &5864 a_lin(i,jsl, njsc(ji)) * mc(ji,jsl,ins) + b_lin(i,jsl,njsc(ji))) ! in mm/d5799 a(ji,jsl) = a_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d 5800 b(ji,jsl) = b_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d 5801 d(ji,jsl) = d_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm^2/d 5802 k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,ji), & 5803 a_lin(i,jsl,ji) * mc(ji,jsl,ins) + b_lin(i,jsl,ji)) ! in mm/d 5865 5804 END DO 5866 5805 END DO … … 5887 5826 !_ hydrol_soil_froz 5888 5827 5889 SUBROUTINE hydrol_soil_froz( kjpindex,ins,njsc)5828 SUBROUTINE hydrol_soil_froz(nvan, avan, mcr, mcs, kjpindex,ins,njsc) 5890 5829 5891 5830 IMPLICIT NONE … … 5894 5833 5895 5834 !! 0.1 Input variables 5835 !salma: added the following soil variables 5836 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: nvan !! Van Genuchten coeficients n (unitless) 5837 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: avan !! Van Genuchten coeficients a (mm-1}) 5838 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 5839 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 5896 5840 5897 5841 INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size … … 5925 5869 DO ji=1,kjpindex 5926 5870 ! Van Genuchten parameter for thermodynamical calculation 5927 m = 1. -1./nvan( njsc(ji))5871 m = 1. -1./nvan(ji) 5928 5872 5929 IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr( njsc(ji))+min_sechiba))) THEN5873 IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr(ji)+min_sechiba))) THEN 5930 5874 ! Linear soil freezing or soil moisture below residual 5931 5875 IF (temp_hydro(ji, jsl).GE.(ZeroCelsius+fr_dT/2.)) THEN … … 5944 5888 (temp_hydro(ji,jsl) .LT. (ZeroCelsius+fr_dT/2.)) ) THEN 5945 5889 ! Factor 2.2 from the PhD of Isabelle Gouttevin 5946 x=MIN(((mcs( njsc(ji))-mcr(njsc(ji))) &5947 *((2.2*1000.*avan( njsc(ji))*(ZeroCelsius+fr_dT/2.-temp_hydro(ji, jsl)) &5948 *lhf/ZeroCelsius/10.)**nvan( njsc(ji))+1.)**(-m)) / &5949 (mc(ji,jsl, ins)-mcr( njsc(ji))),1._r_std)5890 x=MIN(((mcs(ji)-mcr(ji)) & 5891 *((2.2*1000.*avan(ji)*(ZeroCelsius+fr_dT/2.-temp_hydro(ji, jsl)) & 5892 *lhf/ZeroCelsius/10.)**nvan(ji)+1.)**(-m)) / & 5893 (mc(ji,jsl, ins)-mcr(ji)),1._r_std) 5950 5894 ELSE 5951 5895 x=0._r_std … … 5955 5899 profil_froz_hydro_ns(ji,jsl,ins) = 1._r_std-x 5956 5900 5957 mc_ns(ji,jsl)=mc(ji,jsl,ins)/mcs( njsc(ji))5901 mc_ns(ji,jsl)=mc(ji,jsl,ins)/mcs(ji) 5958 5902 5959 5903 ENDDO ! loop on grid … … 6397 6341 !_ hydrol_diag_soil 6398 6342 6399 SUBROUTINE hydrol_diag_soil (k jpindex, veget_max, soiltile, njsc, runoff, drainage, &6343 SUBROUTINE hydrol_diag_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, runoff, drainage, & 6400 6344 & evapot, vevapnu, returnflow, reinfiltration, irrigation, & 6401 6345 & shumdiag,shumdiag_perma, k_litt, litterhumdiag, humrel, vegstress, drysoil_frac, tot_melt) … … 6406 6350 6407 6351 !! 0.1 Input variables 6352 !salma: added the following soil variables 6353 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ks !! Hydraulic conductivity at saturation (mm {-1}) 6354 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: nvan !! Van Genuchten coeficients n (unitless) 6355 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: avan !! Van Genuchten coeficients a (mm-1}) 6356 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcr !! Residual volumetric water content (m^{3} m^{-3}) 6357 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated volumetric water content (m^{3} m^{-3}) 6358 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3}) 6359 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3}) 6408 6360 6409 6361 ! input scalar … … 6558 6510 i= MAX(MIN(INT((imax-imin)*tmc_litter_ratio)+imin, imax-1), imin) 6559 6511 ENDIF 6560 k_tmp = MAX(k_lin(i,1, njsc(ji))*ks(njsc(ji)), zero)6512 k_tmp = MAX(k_lin(i,1,ji)*ks(ji), zero) 6561 6513 k_litt(ji) = k_litt(ji) + vegtot(ji)*soiltile(ji,jst) * SQRT(k_tmp) ! grid-cell average 6562 6514 ENDDO … … 6630 6582 DO ji=1,kjpindex 6631 6583 shumdiag(ji,jsl) = shumdiag(ji,jsl) + soil_wet_ns(ji,jsl,jst) * soiltile(ji,jst) * & 6632 ((mcs( njsc(ji))-mcw(njsc(ji)))/(mcfc(njsc(ji))-mcw(njsc(ji))))6584 ((mcs(ji)-mcw(ji))/(mcfc(ji)-mcw(ji))) 6633 6585 shumdiag(ji,jsl) = MAX(MIN(shumdiag(ji,jsl), un), zero) 6634 6586 ENDDO … … 6641 6593 DO jsl=1,nslm 6642 6594 DO ji=1,kjpindex 6643 shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs( njsc(ji)))6595 shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(ji)) 6644 6596 shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero) 6645 6597 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
