Changeset 7206
- Timestamp:
- 2021-05-31T16:50:43+02:00 (4 years ago)
- Location:
- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba
- Files:
-
- 2 edited
-
sechiba.f90 (modified) (2 diffs)
-
thermosoil.f90 (modified) (29 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90
r6954 r7206 497 497 498 498 !! 1.10 Initialization of soil thermodynamics 499 CALL thermosoil_initialize (kjit, kjpindex, rest_id, &499 CALL thermosoil_initialize (kjit, kjpindex, rest_id, mcs, & 500 500 temp_sol_new, snow, shumdiag_perma, & 501 501 soilcap, soilflx, stempdiag, & … … 771 771 !! 7. Compute soil thermodynamics 772 772 CALL thermosoil_main (kjit, kjpindex, & 773 index, indexgrnd, &773 index, indexgrnd, mcs, & 774 774 temp_sol_new, snow, soilcap, soilflx, & 775 775 shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, & -
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90
r7199 r7206 152 152 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz5 !! Used for numerical calculation [-] 153 153 !$OMP THREADPRIVATE(dz5) 154 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcs !! Saturation humidity [m3/m3]155 !$OMP THREADPRIVATE(mcs)156 154 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: QZ !! quartz content [-] 157 155 !$OMP THREADPRIVATE(QZ) … … 231 229 !! \n 232 230 !_ ============================================================================================================================== 233 SUBROUTINE thermosoil_initialize (kjit, kjpindex, rest_id, 231 SUBROUTINE thermosoil_initialize (kjit, kjpindex, rest_id, mcs, & 234 232 temp_sol_new, snow, shumdiag_perma, & 235 233 soilcap, soilflx, stempdiag, & … … 244 242 INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) 245 243 INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (unitless) 244 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) 246 245 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step, 247 246 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg) … … 370 369 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of tmc_layt','','') 371 370 372 ALLOCATE (mcs(nscm),stat=ier)373 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcs','','')374 375 371 ALLOCATE (QZ(nscm),stat=ier) 376 372 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of QZ','','') … … 379 375 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_dry_ns','','') 380 376 381 ALLOCATE (so_capa_ice( nscm),stat=ier)377 ALLOCATE (so_capa_ice(kjpindex),stat=ier) 382 378 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_ice','','') 383 379 … … 388 384 QZ(:) = QZ_fao(:) 389 385 so_capa_dry_ns(:) = so_capa_dry_ns_fao(:) 390 mcs(:) = mcs_fao(:)391 386 CASE (13) !Salma changed from 12 to 13 for the new class Oxisols 392 387 QZ(:) = QZ_usda(:) 393 388 so_capa_dry_ns(:) = so_capa_dry_ns_usda(:) 394 mcs(:) = mcs_usda(:)395 389 CASE DEFAULT 396 390 WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler, fao and usda according to the map' … … 504 498 505 499 ! Calculate so_capa_ice 506 DO jsc = 1, nscm 507 so_capa_ice(jsc) = so_capa_dry + mcs(jsc)*capa_ice*rho_ice 508 END DO 500 so_capa_ice(:) = so_capa_dry + mcs(:)*capa_ice*rho_ice 509 501 IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice(:)=', so_capa_ice(:),' and capa_ice=',capa_ice 510 502 … … 530 522 CALL thermosoil_coef (& 531 523 kjpindex, temp_sol_new, snow, njsc, & 532 frac_snow_veg, frac_snow_nobio, totfrac_nobio, &524 mcs, frac_snow_veg, frac_snow_nobio, totfrac_nobio, & 533 525 snowdz, snowrho, snowtemp, pb, & 534 526 ptn, & … … 582 574 583 575 SUBROUTINE thermosoil_main (kjit, kjpindex, & 584 index, indexgrnd, &576 index, indexgrnd, mcs, & 585 577 temp_sol_new, snow, soilcap, soilflx, & 586 578 shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, & … … 601 593 INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd !! Indeces of the points on the 3D map (vertical 602 594 !! dimension towards the ground) (unitless) 595 REAL(i_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) 603 596 REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! Surface temperature at the present time-step, 604 597 !! Ts @tex ($K$) @endtex … … 756 749 CALL thermosoil_coef (& 757 750 kjpindex, temp_sol_new, snow, njsc, & 758 frac_snow_veg, frac_snow_nobio, totfrac_nobio, &751 mcs, frac_snow_veg, frac_snow_nobio, totfrac_nobio, & 759 752 snowdz, snowrho, snowtemp, pb, & 760 753 ptn, & … … 918 911 919 912 SUBROUTINE thermosoil_coef (kjpindex, temp_sol_new, snow, njsc, & 920 frac_snow_veg, frac_snow_nobio, totfrac_nobio,&913 mcs, frac_snow_veg, frac_snow_nobio,totfrac_nobio, & 921 914 snowdz, snowrho, snowtemp, pb, & 922 915 ptn, & … … 932 925 REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! snow mass @tex ($Kg$) @endtex 933 926 INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 927 REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) 934 928 REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area 935 929 REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area … … 994 988 ! Computation of the soil thermal properties; snow properties are also accounted for 995 989 IF (ok_freeze_thermix) THEN 996 CALL thermosoil_getdiff( kjpindex, snow, ptn, njsc, snowrho, snowtemp, pb )990 CALL thermosoil_getdiff( kjpindex, snow, ptn, mcs, njsc, snowrho, snowtemp, pb ) 997 991 ELSE 998 992 ! Special case without soil freezing 999 CALL thermosoil_getdiff_old_thermix_without_snow( kjpindex, njsc, snowrho, snowtemp, pb )993 CALL thermosoil_getdiff_old_thermix_without_snow( kjpindex, mcs, njsc, snowrho, snowtemp, pb ) 1000 994 ENDIF 1001 995 … … 1258 1252 !================================================================================================================================ 1259 1253 1260 SUBROUTINE thermosoil_cond (kjpindex, njsc, smc, qz, sh2o, cnd)1254 SUBROUTINE thermosoil_cond (kjpindex, njsc, mcs, smc, qz, sh2o, cnd) 1261 1255 1262 1256 !! 0. Variables and parameter declaration … … 1265 1259 INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) 1266 1260 INTEGER(i_std), DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 1261 REAL(r_std), DIMENSION (kjpindex), INTENT(IN) :: mcs !! Saturated moisture content (m3/m3) 1267 1262 REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3) 1268 1263 REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1) … … 1297 1292 1298 1293 !! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1) 1299 gammd = (1. - mcs(j st))*2700.1294 gammd = (1. - mcs(ji))*2700. 1300 1295 thkdry = (0.135* gammd+ 64.7)/ (2700. - 0.947* gammd) 1301 1296 … … 1312 1307 DO jg = 1,ngrnd 1313 1308 !! 1.4. saturation ratio 1314 satratio(ji,jg) = smc(ji,jg) / mcs(j st)1309 satratio(ji,jg) = smc(ji,jg) / mcs(ji) 1315 1310 1316 1311 !! 1.5. Saturated Thermal Conductivity (thksat) 1317 1312 IF ( smc(ji,jg) > min_sechiba ) THEN 1318 1313 xunfroz(ji,jg) = sh2o(ji,jg) / smc(ji,jg) ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen)) 1319 xu(ji,jg) = xunfroz(ji,jg) * mcs(j st) ! Unfrozen volume for saturation (porosity*xunfroz)1320 thksat(ji,jg) = thks ** (1. - mcs(j st))* THKICE ** (mcs(jst) - xu(ji,jg))* THKW ** (xu(ji,jg))1314 xu(ji,jg) = xunfroz(ji,jg) * mcs(ji) ! Unfrozen volume for saturation (porosity*xunfroz) 1315 thksat(ji,jg) = thks ** (1. - mcs(ji))* THKICE ** (mcs(ji) - xu(ji,jg))* THKW ** (xu(ji,jg)) 1321 1316 ELSE 1322 1317 ! this value will not be used since ake=0 for this case … … 1606 1601 !_ ================================================================================================================================ 1607 1602 1608 SUBROUTINE thermosoil_getdiff( kjpindex, snow, ptn, njsc, snowrho, snowtemp, pb )1603 SUBROUTINE thermosoil_getdiff( kjpindex, snow, ptn, mcs, njsc, snowrho, snowtemp, pb ) 1609 1604 1610 1605 !! 0. Variables and parameter declaration … … 1614 1609 REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow !! Snow mass 1615 1610 INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 1611 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) 1616 1612 REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density 1617 1613 REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K) 1618 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pres ure (hPa)1614 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa) 1619 1615 REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(in) :: ptn !! Soil temperature profile 1620 1616 … … 1640 1636 DO jg = 1, ngrnd 1641 1637 jst = njsc(ji) 1642 pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(j st)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)1638 pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg) 1643 1639 ! 1644 1640 ! 2. Calculate volumetric heat capacity with allowance for permafrost … … 1651 1647 profil_froz(ji,jg) = 1. 1652 1648 pcappa_supp(ji,jg)= 0. 1653 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(j st)) + so_capa_ice(jst) * tmc_layt(ji,jg) / mille / dlt(jg)1654 rho_tot = rho_soil * (1-mcs(j st)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg)1649 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + so_capa_ice(ji) * tmc_layt(ji,jg) / mille / dlt(jg) 1650 rho_tot = rho_soil * (1-mcs(ji)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg) 1655 1651 pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot 1656 1652 … … 1660 1656 profil_froz(ji,jg) = 0. 1661 1657 pcappa_supp(ji,jg)= 0. 1662 rho_tot = rho_soil * (1-mcs(j st)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg)1658 rho_tot = rho_soil * (1-mcs(ji)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg) 1663 1659 pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot 1664 1660 ELSE … … 1668 1664 1669 1665 IF (ok_freeze_thaw_latent_heat) THEN 1670 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(j st)) + &1666 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + & 1671 1667 water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + & 1672 so_capa_ice(j st) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + &1673 shum_ngrnd_perma(ji,jg)*mcs(j st)*lhf*rho_water/fr_dT1668 so_capa_ice(ji) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + & 1669 shum_ngrnd_perma(ji,jg)*mcs(ji)*lhf*rho_water/fr_dT 1674 1670 ELSE 1675 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(j st)) + &1671 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + & 1676 1672 water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + & 1677 so_capa_ice(j st) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)1673 so_capa_ice(ji) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) 1678 1674 ENDIF 1679 1675 1680 rho_tot = rho_soil* (1-mcs(j st)) + &1676 rho_tot = rho_soil* (1-mcs(ji)) + & 1681 1677 rho_water * tmc_layt(ji,jg)/mille / dlt(jg) * xx + & 1682 1678 rho_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) 1683 1679 pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot 1684 1680 1685 pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)*mcs(j st)*lhf*rho_water/fr_dT*zx2(ji,jg)*dlt(jg)1681 pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)*mcs(ji)*lhf*rho_water/fr_dT*zx2(ji,jg)*dlt(jg) 1686 1682 1687 1683 ENDIF … … 1710 1706 ! 1711 1707 IF (ok_freeze_thaw_latent_heat) THEN 1712 CALL thermosoil_cond (kjpindex, njsc, mc _layt, QZ, mcl_layt*(1-profil_froz), pkappa)1708 CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt*(1-profil_froz), pkappa) 1713 1709 ELSE 1714 CALL thermosoil_cond (kjpindex, njsc, mc _layt, QZ, mcl_layt, pkappa)1710 CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt, pkappa) 1715 1711 ENDIF 1716 1712 … … 1744 1740 !_ ================================================================================================================================ 1745 1741 1746 SUBROUTINE thermosoil_getdiff_old_thermix_without_snow( kjpindex, njsc, snowrho, snowtemp, pb )1742 SUBROUTINE thermosoil_getdiff_old_thermix_without_snow( kjpindex, mcs, njsc, snowrho, snowtemp, pb ) 1747 1743 1748 1744 !! 0. Variables and parameter declaration … … 1751 1747 INTEGER(i_std), INTENT(in) :: kjpindex 1752 1748 INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 1749 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) 1753 1750 REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density 1754 1751 REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K) 1755 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pres ure (hPa)1752 REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa) 1756 1753 1757 1754 … … 1764 1761 DO ji = 1,kjpindex 1765 1762 jst = njsc(ji) 1766 pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(j st)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)1763 pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg) 1767 1764 pcapa(ji,jg) = pcapa_tmp(ji, jg) 1768 1765 pcapa_en(ji,jg) = pcapa_tmp(ji, jg) … … 1770 1767 ENDDO 1771 1768 1772 CALL thermosoil_cond (kjpindex, njsc, mc _layt, QZ, mcl_layt, pkappa)1769 CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt, pkappa) 1773 1770 1774 1771 IF (brk_flag == 1) THEN
Note: See TracChangeset
for help on using the changeset viewer.
