Changeset 7199 for branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba
- Timestamp:
- 2021-05-27T14:35:32+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90
r6954 r7199 145 145 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: profil_froz !! Frozen fraction of the soil on hydrological levels (-) 146 146 !$OMP THREADPRIVATE(profil_froz) 147 REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat !! Accumulated latent heat for the whole soil (J) 147 REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat !! Latent heat released or consumed in the freezing/thawing processes summed vertically 148 !! for the whole soil (J/m2) and on the whole simulation to check/correct energy conservation 148 149 !$OMP THREADPRIVATE(e_soil_lat) 149 REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp !! Additional heat capacity due to soil freezing for each soil layer (J/K)150 REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp !! Additional surfacic heat capacity due to soil freezing for each soil layer (J/K/m2) 150 151 !$OMP THREADPRIVATE(pcappa_supp) 151 152 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz5 !! Used for numerical calculation [-] … … 153 154 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcs !! Saturation humidity [m3/m3] 154 155 !$OMP THREADPRIVATE(mcs) 155 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: SMCMAX !! Soil porosity [m3/m3]156 !$OMP THREADPRIVATE(SMCMAX)157 156 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: QZ !! quartz content [-] 158 157 !$OMP THREADPRIVATE(QZ) 159 158 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_dry_ns !! Dry soil Heat capacity of soils,J.m^{-3}.K^{-1} 160 159 !$OMP THREADPRIVATE(so_capa_dry_ns) 160 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_ice !! Heat capacity of saturated frozen soil (J/K/m3) 161 !$OMP THREADPRIVATE(so_capa_ice) 161 162 REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_layt !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface 162 163 !$OMP THREADPRIVATE(mc_layt) … … 271 272 272 273 !! 0.4 Local variables 273 INTEGER(i_std) :: ier, i, jg 274 INTEGER(i_std) :: ier, i, jg, jsc 274 275 LOGICAL :: calculate_coef !! Local flag to initialize variables by call to thermosoil_coef 275 276 !_ ================================================================================================================================ … … 372 373 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcs','','') 373 374 374 ALLOCATE (SMCMAX(nscm),stat=ier)375 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of SMCMAX','','')376 377 375 ALLOCATE (QZ(nscm),stat=ier) 378 376 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of QZ','','') … … 380 378 ALLOCATE (so_capa_dry_ns(nscm),stat=ier) 381 379 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_dry_ns','','') 380 381 ALLOCATE (so_capa_ice(nscm),stat=ier) 382 IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_ice','','') 382 383 383 384 … … 385 386 SELECTCASE (nscm) 386 387 CASE (3) 387 SMCMAX(:) = SMCMAX_fao(:)388 388 QZ(:) = QZ_fao(:) 389 389 so_capa_dry_ns(:) = so_capa_dry_ns_fao(:) 390 390 mcs(:) = mcs_fao(:) 391 391 CASE (13) !Salma changed from 12 to 13 for the new class Oxisols 392 SMCMAX(:) = SMCMAX_usda(:) 393 QZ(:) = QZ_usda(:) 392 QZ(:) = QZ_usda(:) 394 393 so_capa_dry_ns(:) = so_capa_dry_ns_usda(:) 395 394 mcs(:) = mcs_usda(:) … … 505 504 506 505 ! Calculate so_capa_ice 507 so_capa_ice = so_capa_dry + poros*capa_ice*rho_ice 508 IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice=', so_capa_ice,' using poros=',poros,' and capa_ice=',capa_ice 506 DO jsc = 1, nscm 507 so_capa_ice(jsc) = so_capa_dry + mcs(jsc)*capa_ice*rho_ice 508 END DO 509 IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice(:)=', so_capa_ice(:),' and capa_ice=',capa_ice 509 510 510 511 ! Computing some usefull constants for the numerical scheme … … 1257 1258 !================================================================================================================================ 1258 1259 1259 SUBROUTINE thermosoil_cond (kjpindex, njsc, smc, qz, s mcmax, sh2o, cnd)1260 SUBROUTINE thermosoil_cond (kjpindex, njsc, smc, qz, sh2o, cnd) 1260 1261 1261 1262 !! 0. Variables and parameter declaration … … 1265 1266 INTEGER(i_std), DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 1266 1267 REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3) 1267 REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1) 1268 REAL(r_std), DIMENSION (nscm), INTENT(IN) :: smcmax !! Soil Porosity (0-1) 1268 REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1) 1269 1269 REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o 1270 1270 … … 1297 1297 1298 1298 !! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1) 1299 gammd = (1. - smcmax(jst))*2700.1299 gammd = (1. - mcs(jst))*2700. 1300 1300 thkdry = (0.135* gammd+ 64.7)/ (2700. - 0.947* gammd) 1301 1301 … … 1312 1312 DO jg = 1,ngrnd 1313 1313 !! 1.4. saturation ratio 1314 satratio(ji,jg) = smc(ji,jg) / smcmax(jst)1314 satratio(ji,jg) = smc(ji,jg) / mcs(jst) 1315 1315 1316 1316 !! 1.5. Saturated Thermal Conductivity (thksat) 1317 1317 IF ( smc(ji,jg) > min_sechiba ) THEN 1318 1318 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) * smcmax(jst) ! Unfrozen volume for saturation (porosity*xunfroz)1320 thksat(ji,jg) = thks ** (1. - smcmax(jst))* THKICE ** (smcmax(jst) - xu(ji,jg))* THKW ** (xu(ji,jg))1319 xu(ji,jg) = xunfroz(ji,jg) * mcs(jst) ! Unfrozen volume for saturation (porosity*xunfroz) 1320 thksat(ji,jg) = thks ** (1. - mcs(jst))* THKICE ** (mcs(jst) - xu(ji,jg))* THKW ** (xu(ji,jg)) 1321 1321 ELSE 1322 1322 ! this value will not be used since ake=0 for this case … … 1622 1622 REAL :: xx !! Unfrozen fraction of the soil 1623 1623 REAL(r_std), DIMENSION(kjpindex) :: snow_h 1624 REAL(r_std), DIMENSION(kjpindex,ngrnd) :: zx1, zx2 1625 REAL :: cap_iw !! Heat capacity of ice/water mixture 1624 REAL(r_std), DIMENSION(kjpindex,ngrnd) :: zx1, zx2 1626 1625 INTEGER :: ji,jg 1627 1626 INTEGER :: jst … … 1652 1651 profil_froz(ji,jg) = 1. 1653 1652 pcappa_supp(ji,jg)= 0. 1654 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + so_capa_ice * tmc_layt(ji,jg) / mille / dlt(jg)1653 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + so_capa_ice(jst) * tmc_layt(ji,jg) / mille / dlt(jg) 1655 1654 rho_tot = rho_soil * (1-mcs(jst)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg) 1656 1655 pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot … … 1668 1667 profil_froz(ji,jg) = (1. - xx) 1669 1668 1670 ! net heat capacity of the ice/water mixture1671 cap_iw = xx * so_capa_wet + (1.-xx) * so_capa_ice1672 1669 IF (ok_freeze_thaw_latent_heat) THEN 1673 1670 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + & 1674 1671 water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + & 1675 so_capa_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + &1676 shum_ngrnd_perma(ji,jg)* poros*lhf*rho_water/fr_dT1672 so_capa_ice(jst) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + & 1673 shum_ngrnd_perma(ji,jg)*mcs(jst)*lhf*rho_water/fr_dT 1677 1674 ELSE 1678 1675 pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + & 1679 1676 water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + & 1680 so_capa_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)1677 so_capa_ice(jst) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) 1681 1678 ENDIF 1682 1679 … … 1686 1683 pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot 1687 1684 1688 pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)* poros*lhf*rho_water/fr_dT*zx2(ji,jg)*dlt(jg)1685 pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)*mcs(jst)*lhf*rho_water/fr_dT*zx2(ji,jg)*dlt(jg) 1689 1686 1690 1687 ENDIF … … 1713 1710 ! 1714 1711 IF (ok_freeze_thaw_latent_heat) THEN 1715 CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, SMCMAX,mcl_layt*(1-profil_froz), pkappa)1712 CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt*(1-profil_froz), pkappa) 1716 1713 ELSE 1717 CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, SMCMAX,mcl_layt, pkappa)1714 CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt, pkappa) 1718 1715 ENDIF 1719 1716 … … 1773 1770 ENDDO 1774 1771 1775 CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, SMCMAX,mcl_layt, pkappa)1772 CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt, pkappa) 1776 1773 1777 1774 IF (brk_flag == 1) THEN
Note: See TracChangeset
for help on using the changeset viewer.
