SUBROUTINE SOIL_HTC ( 2,6SOILHT7A.23
& NPNTS,NSHYD,SOIL_PTS,SOIL_INDEX,B,DZ,HCAP,HCON,SATHH SOILHT7A.24
&,SNOW_SOIL_HTF,SOIL_SURF_HTF,TIMESTEP,V_SAT,W_FLUX ARE1F405.43
&,SMCL,STHU,STHF,TSOIL SOILHT7A.26
+,LTIMER SOILHT7A.27
+) SOILHT7A.28
SOILHT7A.29
IMPLICIT NONE SOILHT7A.30
! SOILHT7A.31
! Description: SOILHT7A.32
! Updates deep soil temperatures, frozen and unfrozen SOILHT7A.33
! frozen soil water content. Calls the following subroutines: SOILHT7A.34
! SOILHT7A.35
! HEAT_CON - to calculate the soil thermal conductivity SOILHT7A.36
! (Cox, 6/95) SOILHT7A.37
! SOILHT7A.38
! Documentation : UM Documentation Paper 25 SOILHT7A.39
! SOILHT7A.40
! Current Code Owner : David Gregory SOILHT7A.41
! SOILHT7A.42
! History: SOILHT7A.43
! Version Date Comment SOILHT7A.44
! ------- ---- ------- SOILHT7A.45
! 4.1 New deck. Peter Cox SOILHT7A.46
! 4.4 7/97 MOSES II. Richard Essery SOILHT7A.47
! 4.5 26/5/98 Correction to stop overflow when SMCL very small. ACB2F405.13
! C.Bunton ACB2F405.14
!LL 4.5 18/06/98 Changed Timer calls to indicate non-barrier GPB8F405.39
!LL P.Burton GPB8F405.40
! SOILHT7A.48
! Code Description: SOILHT7A.49
! Language: FORTRAN 77 + common extensions. SOILHT7A.50
! SOILHT7A.51
! System component covered: P25 SOILHT7A.52
! System Task: P25 SOILHT7A.53
! SOILHT7A.54
SOILHT7A.55
! Global variables: SOILHT7A.56
*CALL C_SOILH 
SOILHT7A.57
SOILHT7A.58
! Subroutine arguments SOILHT7A.59
! Scalar arguments with intent(IN) : SOILHT7A.60
INTEGER SOILHT7A.61
& NPNTS ! IN Number of gridpoints. SOILHT7A.62
&,NSHYD ! IN Number of soil moisture levels. SOILHT7A.63
&,SOIL_PTS ! IN Number of soil points. SOILHT7A.64
SOILHT7A.65
REAL SOILHT7A.66
& TIMESTEP ! IN Model timestep (s). SOILHT7A.67
SOILHT7A.68
SOILHT7A.69
! Array arguments with intent(IN) : SOILHT7A.70
INTEGER SOILHT7A.71
& SOIL_INDEX(NPNTS) ! IN Array of soil points. SOILHT7A.72
SOILHT7A.73
REAL SOILHT7A.74
& B(NPNTS) ! IN Clapp-Hornberger exponent. SOILHT7A.75
&,DZ(NSHYD) ! IN Thicknesses of the soil layers (m). SOILHT7A.76
&,HCAP(NPNTS) ! IN Soil heat capacity (J/K/m3). SOILHT7A.77
&,HCON(NPNTS) ! IN Soil thermal conductivity (W/m/K). SOILHT7A.78
&,SATHH(NPNTS) ! IN Saturated soil water pressure (m). SOILHT7A.79
&,SMCL(NPNTS,NSHYD) ! IN Soil moisture content of each SOILHT7A.80
! ! layer (kg/m2). SOILHT7A.81
&,SNOW_SOIL_HTF(NPNTS) ! IN Heat flux from snow to soil (W/m2). SOILHT7A.83
&,SOIL_SURF_HTF(NPNTS) ! IN Net downward surface heat flux (W/m2). SOILHT7A.84
&,V_SAT(NPNTS) ! IN Volumetric soil moisture SOILHT7A.85
! ! concentration at saturation SOILHT7A.86
! ! (m3 H2O/m3 soil). SOILHT7A.87
+,W_FLUX(NPNTS,0:NSHYD)! IN The fluxes of water between layers SOILHT7A.88
! ! (kg/m2/s). SOILHT7A.89
C SOILHT7A.90
LOGICAL LTIMER ! Logical switch for TIMER diags SOILHT7A.91
SOILHT7A.92
! Array arguments with intent(OUT) : SOILHT7A.93
SOILHT7A.94
! Array arguments with intent(INOUT) : SOILHT7A.95
REAL SOILHT7A.96
& STHF(NPNTS,NSHYD) ! INOUT Frozen soil moisture content of SOILHT7A.97
! ! each layer as a fraction of SOILHT7A.98
! ! saturation. SOILHT7A.99
&,STHU(NPNTS,NSHYD) ! INOUT Unfrozen soil moisture content of SOILHT7A.100
! ! each layer as a fraction of SOILHT7A.101
! ! saturation. SOILHT7A.102
&,TSOIL(NPNTS,NSHYD) ! INOUT Sub-surface temperatures (K). SOILHT7A.103
SOILHT7A.104
! Local scalars: SOILHT7A.105
INTEGER SOILHT7A.106
& I,J,M,N ! WORK Loop counters. SOILHT7A.107
&,J_ITER ! WORK Number of soil points which require SOILHT7A.108
! ! iteration. SOILHT7A.109
&,MMAX ! WORK Maximum number of iterations on SOILHT7A.110
! ! temperature. SOILHT7A.111
&,FIRST_SOIL_PT ! First soil point SOILHT7A.112
SOILHT7A.113
REAL SOILHT7A.114
& FACUR ! WORK Required flux conservation accuracy SOILHT7A.115
! ! (W/m2). SOILHT7A.116
&,TACUR ! WORK Required accuracy of temperature SOILHT7A.117
! ! calculation (Celsius). SOILHT7A.118
PARAMETER( MMAX=3, FACUR=0.01, TACUR=0.00000 ) SOILHT7A.119
SOILHT7A.120
! Local arrays: SOILHT7A.121
INTEGER SOILHT7A.122
& ITER_PTS(NPNTS) ! WORK Array of soil points which require SOILHT7A.123
! ! iteration. SOILHT7A.124
SOILHT7A.125
REAL SOILHT7A.126
& CEACUR(NPNTS) ! WORK Flux conservation accuracy of the SOILHT7A.127
! ! calculation (W/m2) SOILHT7A.128
&,DHSL0(NPNTS,NSHYD) ! WORK Total heat increment to the layer SOILHT7A.129
! ! (J/m2/timestep) SOILHT7A.130
&,DHSL(NPNTS,NSHYD) ! WORK The heat available to update the SOILHT7A.131
! ! layer temperature (J/m2/timestep) SOILHT7A.132
&,DHSLA(NPNTS,NSHYD) ! WORK The heat increment due to advection o SOILHT7A.133
! ! heat by moisture fluxes (J/m2/timeste SOILHT7A.134
&,DSMCLF(NPNTS,NSHYD) ! WORK The increment to the layer frozen SOILHT7A.135
! ! soil moisture (kg/m2/timestep). SOILHT7A.136
&,DTHU(NPNTS,NSHYD) ! WORK Rate of change of volumetric unfrozen SOILHT7A.137
! ! soil moisture concentration with SOILHT7A.138
! ! temperature (m3 liquid H2O/m3 soil/K) SOILHT7A.139
&,DTSL(NPNTS,NSHYD) ! WORK The increment to the layer temperatur SOILHT7A.140
! ! (K/timestep). SOILHT7A.141
&,HCAPT(NPNTS) ! WORK The total volumetric heat capacity SOILHT7A.142
! ! (soil+water) of the layer (J/m3/K). SOILHT7A.143
&,HC(NPNTS,NSHYD) ! WORK The thermal conductivity of each SOILHT7A.144
! ! layer (W/m/K). SOILHT7A.145
&,HCONS(NPNTS) ! WORK The thermal conductivity between SOILHT7A.146
! ! adjacent soil layers (W/m/K). SOILHT7A.147
&,H_FLUX(NPNTS,0:NSHYD) !WORK The fluxes of heat between layers SOILHT7A.148
! ! (W/m2). SOILHT7A.149
&,SMCLF(NPNTS,NSHYD) ! WORK Frozen moisture content of each SOILHT7A.150
! ! soil layer (kg/m2). SOILHT7A.151
&,SMCLF0(NPNTS,NSHYD) ! WORK Previous value of SMCLF (kg/m2). SOILHT7A.152
&,SMCLSAT(NPNTS,NSHYD) ! WORK The saturation moisture content of SOILHT7A.153
! ! each layer (kg/m2). SOILHT7A.154
&,SMCLU(NPNTS,NSHYD) ! WORK Unfrozen moisture content of each SOILHT7A.155
! ! soil layer (kg/m2). SOILHT7A.156
&,SMCLU0(NPNTS,NSHYD) ! WORK Previous value of SMCLU (kg/m2). SOILHT7A.157
&,SMCFU(NPNTS) ! WORK Fractional saturation (unfrozen water SOILHT7A.158
! ! at layer boundaries. SOILHT7A.159
&,SMCFF(NPNTS) ! WORK Fractional saturation (frozen water) SOILHT7A.160
! ! at layer boundaries. SOILHT7A.161
&,TMAX(NPNTS,NSHYD) ! WORK Temperature above which all water is SOILHT7A.162
! ! unfrozen (Celsius) SOILHT7A.163
&,TSL(NPNTS,0:NSHYD) ! WORK Soil layer temperatures (Celsius) SOILHT7A.164
! ! TSL(0) temperature of incoming water. SOILHT7A.165
! ! TSL(1:NSHYD) sub-surface soil SOILHT7A.166
! ! temperatures . SOILHT7A.167
&,TSL0(NPNTS,0:NSHYD) ! WORK Previous value of TSL (Celsius). SOILHT7A.168
SOILHT7A.169
! Function & Subroutine calls: SOILHT7A.170
EXTERNAL SOILHT7A.171
& HEAT_CON SOILHT7A.172
SOILHT7A.173
*CALL C_DENSTY SOILHT7A.174
*CALL C_LHEAT SOILHT7A.175
*CALL C_PERMA SOILHT7A.176
*CALL C_0_DG_C SOILHT7A.177
SOILHT7A.178
IF (LTIMER) THEN SOILHT7A.179
CALL TIMER('SOILHTC ',103) GPB8F405.41
ENDIF SOILHT7A.181
SOILHT7A.182
!----------------------------------------------------------------------- SOILHT7A.183
! Initialisations SOILHT7A.184
!----------------------------------------------------------------------- SOILHT7A.185
FIRST_SOIL_PT=SOIL_INDEX(1) SOILHT7A.186
SOILHT7A.187
DO J=1,SOIL_PTS SOILHT7A.188
I=SOIL_INDEX(J) SOILHT7A.189
TSL(I,0)=TSOIL(I,1)-ZERODEGC SOILHT7A.190
ENDDO SOILHT7A.191
SOILHT7A.192
DO N=1,NSHYD SOILHT7A.193
DO J=1,SOIL_PTS SOILHT7A.194
I=SOIL_INDEX(J) SOILHT7A.195
!----------------------------------------------------------------------- SOILHT7A.196
! Define soil layer temperatures TSL (in celsius). SOILHT7A.197
!----------------------------------------------------------------------- SOILHT7A.198
TSL(I,N)=TSOIL(I,N)-ZERODEGC SOILHT7A.199
ENDDO SOILHT7A.200
ENDDO SOILHT7A.201
SOILHT7A.202
DO N=1,NSHYD SOILHT7A.203
! CDIR$ IVDEP here would force vectorization but changes results! SOILHT7A.204
DO J=1,SOIL_PTS SOILHT7A.205
I=SOIL_INDEX(J) SOILHT7A.206
!----------------------------------------------------------------------- SOILHT7A.207
! Diagnose the frozen and unfrozen water. SOILHT7A.208
!----------------------------------------------------------------------- SOILHT7A.209
SOILHT7A.210
SMCLSAT(I,N)=RHO_WATER*DZ(N)*V_SAT(I) SOILHT7A.211
SMCLF(I,N)=SMCLSAT(I,N)*STHF(I,N) SOILHT7A.212
SMCLU(I,N)=SMCL(I,N)-SMCLF(I,N) SOILHT7A.213
ENDDO ! J=1,SOIL_PTS SOILHT7A.214
ENDDO ! N=1,NSHYD SOILHT7A.215
SOILHT7A.216
!---------------------------------------------------------------------- SOILHT7A.217
! Calculate the heat conductivity between adjacent layers SOILHT7A.218
!---------------------------------------------------------------------- SOILHT7A.219
SOILHT7A.220
DO N=1,NSHYD-1 SOILHT7A.221
DO J=1,SOIL_PTS SOILHT7A.222
I=SOIL_INDEX(J) SOILHT7A.223
SMCFU(I)=(DZ(N+1)*STHU(I,N)+DZ(N)*STHU(I,N+1)) SOILHT7A.224
& /(DZ(N+1)+DZ(N)) SOILHT7A.225
SMCFF(I)=(DZ(N+1)*STHF(I,N)+DZ(N)*STHF(I,N+1)) SOILHT7A.226
& /(DZ(N+1)+DZ(N)) SOILHT7A.227
ENDDO SOILHT7A.228
CALL HEAT_CON(NPNTS-FIRST_SOIL_PT+1,HCON(FIRST_SOIL_PT) SOILHT7A.229
&, SMCFU(FIRST_SOIL_PT),SMCFF(FIRST_SOIL_PT) GJC0F405.66
&, V_SAT(FIRST_SOIL_PT),HCONS(FIRST_SOIL_PT),LTIMER) SOILHT7A.231
DO J=1,SOIL_PTS SOILHT7A.232
I=SOIL_INDEX(J) SOILHT7A.233
HC(I,N)=HCONS(I) SOILHT7A.234
ENDDO SOILHT7A.235
ENDDO SOILHT7A.236
SOILHT7A.237
!-------------------------------------------------------------------- SOILHT7A.238
! Calculate heat fluxes across layer boundaries SOILHT7A.239
!-------------------------------------------------------------------- SOILHT7A.240
DO N=1,NSHYD-1 SOILHT7A.241
DO J=1,SOIL_PTS SOILHT7A.242
I=SOIL_INDEX(J) SOILHT7A.243
H_FLUX(I,N)=-HC(I,N)*2.0*(TSL(I,N+1)-TSL(I,N)) SOILHT7A.244
& /(DZ(N+1)+DZ(N)) SOILHT7A.245
ENDDO SOILHT7A.246
ENDDO SOILHT7A.247
CDIR$ IVDEP SOILHT7A.248
DO J=1,SOIL_PTS SOILHT7A.249
I=SOIL_INDEX(J) SOILHT7A.250
H_FLUX(I,NSHYD)=0.0 SOILHT7A.251
H_FLUX(I,0) = SOIL_SURF_HTF(I) + SNOW_SOIL_HTF(I) ARE1F405.44
ENDDO SOILHT7A.254
SOILHT7A.255
!-------------------------------------------------------------------- SOILHT7A.256
! Calculate the advection of heat by moisture fluxes SOILHT7A.257
!-------------------------------------------------------------------- SOILHT7A.258
DO N=2,NSHYD-1 SOILHT7A.259
DO J=1,SOIL_PTS SOILHT7A.260
I=SOIL_INDEX(J) SOILHT7A.261
DHSLA(I,N)=TIMESTEP*HCAPW*DZ(N)* SOILHT7A.262
& (W_FLUX(I,N-1)*(TSL(I,N-1)-TSL(I,N))/(DZ(N)+DZ(N-1)) SOILHT7A.263
& +W_FLUX(I,N)*(TSL(I,N)-TSL(I,N+1))/(DZ(N)+DZ(N+1))) SOILHT7A.264
ENDDO SOILHT7A.265
ENDDO SOILHT7A.266
SOILHT7A.267
CDIR$ IVDEP SOILHT7A.268
DO J=1,SOIL_PTS SOILHT7A.269
I=SOIL_INDEX(J) SOILHT7A.270
DHSLA(I,1)=TIMESTEP*HCAPW*DZ(1)* SOILHT7A.271
& (W_FLUX(I,0)*(TSL(I,0)-TSL(I,1))/DZ(1) SOILHT7A.272
& +W_FLUX(I,1)*(TSL(I,1)-TSL(I,2))/(DZ(1)+DZ(2))) SOILHT7A.273
DHSLA(I,NSHYD)=TIMESTEP*HCAPW*DZ(NSHYD)* SOILHT7A.274
& W_FLUX(I,NSHYD-1)*(TSL(I,NSHYD-1)-TSL(I,NSHYD)) SOILHT7A.275
& /(DZ(NSHYD)+DZ(NSHYD-1)) SOILHT7A.276
ENDDO SOILHT7A.277
SOILHT7A.278
DO N=1,NSHYD SOILHT7A.279
CDIR$ IVDEP SOILHT7A.280
DO J=1,SOIL_PTS SOILHT7A.281
I=SOIL_INDEX(J) SOILHT7A.282
!----------------------------------------------------------------------- SOILHT7A.283
! Calculate TMAX, the temperature above which all soil water is SOILHT7A.284
! unfrozen. Check that (SMCLSAT/SMCL)**B will not overflow when SMCL is ACB2F405.15
! very small. The function EPSILON gives the number of type (real) of ACB2F405.16
! SMCL that is negligeable compared to 1. ACB2F405.17
!----------------------------------------------------------------------- ACB2F405.18
*IF DEF,SCMA,AND,-DEF,T3E ACB2F405.19
IF (SMCL(I,N) .GT. SMCLSAT(I,N)*1E-20) THEN ACB2F405.20
*ELSE ACB2F405.21
IF (SMCL(I,N) .GT. SMCLSAT(I,N)*(EPSILON(SMCL(I,N)) ) )THEN ACB2F405.22
*ENDIF ACB2F405.23
ACB2F405.24
TMAX(I,N)=-SATHH(I)/DPSIDT SOILHT7A.288
& *(SMCLSAT(I,N)/SMCL(I,N))**(B(I)) SOILHT7A.289
ELSE SOILHT7A.290
TMAX(I,N)=-ZERODEGC SOILHT7A.291
ENDIF SOILHT7A.292
TMAX(I,N)=MAX(TMAX(I,N),-ZERODEGC) SOILHT7A.293
SOILHT7A.294
DHSL0(I,N)=TIMESTEP*(H_FLUX(I,N-1)-H_FLUX(I,N)) SOILHT7A.295
& +DHSLA(I,N) SOILHT7A.296
SOILHT7A.297
SOILHT7A.298
DHSL(I,N)=DHSL0(I,N) SOILHT7A.299
SOILHT7A.300
!----------------------------------------------------------------------- SOILHT7A.301
! Calculate the soil temperature increments SOILHT7A.302
!----------------------------------------------------------------------- SOILHT7A.303
TSL0(I,N)=TSL(I,N) SOILHT7A.304
SMCLF0(I,N)=SMCLF(I,N) SOILHT7A.305
SMCLU0(I,N)=SMCLU(I,N) SOILHT7A.306
SOILHT7A.307
IF (TSL(I,N).GT.TMAX(I,N)) THEN ! All water unfrozen SOILHT7A.308
DTHU(I,N)=0.0 SOILHT7A.309
ELSEIF (TSL(I,N).EQ.TMAX(I,N).AND. ! Remains unfrozen SOILHT7A.310
& DHSL(I,N).GE.0.0) THEN SOILHT7A.311
DTHU(I,N)=0.0 SOILHT7A.312
ELSE ! Phase changes SOILHT7A.313
DTHU(I,N)=DPSIDT*SMCLSAT(I,N) SOILHT7A.314
& /(B(I)*SATHH(I)*RHO_WATER*DZ(N)) SOILHT7A.315
& *(-DPSIDT*TSL(I,N)/SATHH(I))**(-1.0/B(I)-1.0) SOILHT7A.316
ENDIF SOILHT7A.317
SOILHT7A.318
HCAPT(I)=HCAP(I)+(HCAPW-HCAPI)*SMCLU(I,N)/DZ(N) SOILHT7A.319
& +HCAPI*SMCL(I,N)/DZ(N) SOILHT7A.320
& +RHO_WATER*DTHU(I,N)*((HCAPW-HCAPI)*TSL(I,N)+LF) SOILHT7A.321
SOILHT7A.322
DTSL(I,N)=1.0/(HCAPT(I)*DZ(N))*DHSL(I,N) SOILHT7A.323
TSL(I,N)=TSL(I,N)+DTSL(I,N) SOILHT7A.324
SOILHT7A.325
!----------------------------------------------------------------------- SOILHT7A.326
! If the temperature increment is small and frozen water exists SOILHT7A.327
! assume that the excess energy goes into phase change SOILHT7A.328
!----------------------------------------------------------------------- SOILHT7A.329
IF (ABS(DTSL(I,N)).LT.TACUR.AND. SOILHT7A.330
& TSL(I,N).LE.TMAX(I,N)) THEN SOILHT7A.331
DSMCLF(I,N)=-DHSL(I,N)/((HCAPW-HCAPI)*TSL(I,N)+LF) SOILHT7A.332
DSMCLF(I,N)=MAX(DSMCLF(I,N),-SMCLF(I,N)) SOILHT7A.333
DSMCLF(I,N)=MIN(DSMCLF(I,N),SMCLU(I,N)) SOILHT7A.334
SMCLU(I,N)=SMCLU(I,N)-DSMCLF(I,N) SOILHT7A.335
SMCLF(I,N)=SMCLF(I,N)+DSMCLF(I,N) SOILHT7A.336
ENDIF SOILHT7A.337
SOILHT7A.338
!----------------------------------------------------------------------- SOILHT7A.339
! Check to see if the discontinuity in HCAPT at TMAX has been crossed SOILHT7A.340
! - if it has return to TMAX SOILHT7A.341
!----------------------------------------------------------------------- SOILHT7A.342
IF ((TSL(I,N)-TMAX(I,N))*(TSL0(I,N)-TMAX(I,N)) SOILHT7A.343
& .LT.0.0) THEN SOILHT7A.344
DTSL(I,N)=DTSL(I,N)+(TMAX(I,N)-TSL(I,N)) SOILHT7A.345
TSL(I,N)=TMAX(I,N) SOILHT7A.346
ENDIF SOILHT7A.347
SOILHT7A.348
!----------------------------------------------------------------------- SOILHT7A.349
! Diagnose unfrozen and frozen water contents SOILHT7A.350
!----------------------------------------------------------------------- SOILHT7A.351
IF (TSL(I,N).GE.TMAX(I,N)) THEN SOILHT7A.352
SMCLU(I,N)=SMCL(I,N) SOILHT7A.353
SMCLF(I,N)=0.0 SOILHT7A.354
ELSE SOILHT7A.355
SMCLU(I,N)=SMCLSAT(I,N) SOILHT7A.356
& *(-DPSIDT*TSL(I,N)/SATHH(I))**(-1.0/B(I)) SOILHT7A.357
SMCLF(I,N)=SMCL(I,N)-SMCLU(I,N) SOILHT7A.358
ENDIF SOILHT7A.359
SOILHT7A.360
!----------------------------------------------------------------------- SOILHT7A.361
! Calculate the error in heat conservation SOILHT7A.362
!----------------------------------------------------------------------- SOILHT7A.363
DSMCLF(I,N)=SMCLF(I,N)-SMCLF0(I,N) SOILHT7A.364
DHSL(I,N)=DHSL(I,N)-(HCAP(I)*DZ(N)+HCAPW*SMCLU0(I,N) SOILHT7A.365
& +HCAPI*SMCLF0(I,N))*DTSL(I,N) SOILHT7A.366
& -DSMCLF(I,N)*((HCAPI-HCAPW)*TSL0(I,N)-LF) SOILHT7A.367
SOILHT7A.368
!----------------------------------------------------------------------- SOILHT7A.369
! Calculate the error in flux conservation SOILHT7A.370
!----------------------------------------------------------------------- SOILHT7A.371
CEACUR(I)=ABS(DHSL(I,N))/TIMESTEP SOILHT7A.372
SOILHT7A.373
ENDDO SOILHT7A.374
SOILHT7A.375
!----------------------------------------------------------------------- SOILHT7A.376
! Define the array of points which fail to meet the flux criterion SOILHT7A.377
!----------------------------------------------------------------------- SOILHT7A.378
J_ITER=0 SOILHT7A.379
C CDIR$ IVDEP SOILHT7A.380
DO J=1,SOIL_PTS SOILHT7A.381
I=SOIL_INDEX(J) SOILHT7A.382
SOILHT7A.383
IF (CEACUR(I) .GT. FACUR) THEN SOILHT7A.384
J_ITER=J_ITER+1 SOILHT7A.385
ITER_PTS(J_ITER)=J SOILHT7A.386
ENDIF SOILHT7A.387
SOILHT7A.388
ENDDO SOILHT7A.389
SOILHT7A.390
SOILHT7A.391
!----------------------------------------------------------------------- SOILHT7A.392
! Iterate on these points SOILHT7A.393
!----------------------------------------------------------------------- SOILHT7A.394
DO M=1,MMAX-1 SOILHT7A.395
CDIR$ IVDEP SOILHT7A.396
DO J=1,J_ITER SOILHT7A.397
I=SOIL_INDEX(ITER_PTS(J)) SOILHT7A.398
SOILHT7A.399
TSL0(I,N)=TSL(I,N) SOILHT7A.400
SMCLF0(I,N)=SMCLF(I,N) SOILHT7A.401
SMCLU0(I,N)=SMCLU(I,N) SOILHT7A.402
SOILHT7A.403
IF (TSL(I,N).GT.TMAX(I,N)) THEN ! All water unfrozen SOILHT7A.404
DTHU(I,N)=0.0 SOILHT7A.405
ELSEIF (TSL(I,N).EQ.TMAX(I,N).AND. ! Remains unfrozen SOILHT7A.406
& DHSL(I,N).GE.0.0) THEN SOILHT7A.407
DTHU(I,N)=0.0 SOILHT7A.408
ELSE ! Phase changes SOILHT7A.409
DTHU(I,N)=DPSIDT*SMCLSAT(I,N) SOILHT7A.410
& /(B(I)*SATHH(I)*RHO_WATER*DZ(N)) SOILHT7A.411
& *(-DPSIDT*TSL(I,N)/SATHH(I))**(-1.0/B(I)-1.0) SOILHT7A.412
ENDIF SOILHT7A.413
SOILHT7A.414
HCAPT(I)=HCAP(I)+(HCAPW-HCAPI)*SMCLU(I,N)/DZ(N) SOILHT7A.415
& +HCAPI*SMCL(I,N)/DZ(N) SOILHT7A.416
& +RHO_WATER*DTHU(I,N)*((HCAPW-HCAPI)*TSL(I,N)+LF) SOILHT7A.417
SOILHT7A.418
DTSL(I,N)=1.0/(HCAPT(I)*DZ(N))*DHSL(I,N) SOILHT7A.419
TSL(I,N)=TSL(I,N)+DTSL(I,N) SOILHT7A.420
SOILHT7A.421
!----------------------------------------------------------------------- SOILHT7A.422
! If the temperature increment is small and frozen water exists SOILHT7A.423
! assume that the excess energy goes into phase change SOILHT7A.424
!----------------------------------------------------------------------- SOILHT7A.425
IF (ABS(DTSL(I,N)).LT.TACUR.AND. SOILHT7A.426
& TSL(I,N).LE.TMAX(I,N)) THEN SOILHT7A.427
DSMCLF(I,N)=-DHSL(I,N)/((HCAPW-HCAPI)*TSL(I,N)+LF) SOILHT7A.428
DSMCLF(I,N)=MAX(DSMCLF(I,N),-SMCLF(I,N)) SOILHT7A.429
DSMCLF(I,N)=MIN(DSMCLF(I,N),SMCLU(I,N)) SOILHT7A.430
SMCLU(I,N)=SMCLU(I,N)-DSMCLF(I,N) SOILHT7A.431
SMCLF(I,N)=SMCLF(I,N)+DSMCLF(I,N) SOILHT7A.432
ENDIF SOILHT7A.433
SOILHT7A.434
!----------------------------------------------------------------------- SOILHT7A.435
! Check to see if the discontinuity in HCAPT at TMAX has been crossed SOILHT7A.436
! - if it has return to TMAX SOILHT7A.437
!----------------------------------------------------------------------- SOILHT7A.438
IF ((TSL(I,N)-TMAX(I,N))*(TSL0(I,N)-TMAX(I,N)) SOILHT7A.439
& .LT.0.0) THEN SOILHT7A.440
DTSL(I,N)=DTSL(I,N)+(TMAX(I,N)-TSL(I,N)) SOILHT7A.441
TSL(I,N)=TMAX(I,N) SOILHT7A.442
ENDIF SOILHT7A.443
SOILHT7A.444
!----------------------------------------------------------------------- SOILHT7A.445
! Diagnose unfrozen and frozen water contents SOILHT7A.446
!----------------------------------------------------------------------- SOILHT7A.447
IF (TSL(I,N).GE.TMAX(I,N)) THEN SOILHT7A.448
SMCLU(I,N)=SMCL(I,N) SOILHT7A.449
SMCLF(I,N)=0.0 SOILHT7A.450
ELSE SOILHT7A.451
SMCLU(I,N)=SMCLSAT(I,N) SOILHT7A.452
& *(-DPSIDT*TSL(I,N)/SATHH(I))**(-1.0/B(I)) SOILHT7A.453
SMCLF(I,N)=SMCL(I,N)-SMCLU(I,N) SOILHT7A.454
ENDIF SOILHT7A.455
SOILHT7A.456
!----------------------------------------------------------------------- SOILHT7A.457
! Calculate the error in heat conservation SOILHT7A.458
!----------------------------------------------------------------------- SOILHT7A.459
DSMCLF(I,N)=SMCLF(I,N)-SMCLF0(I,N) SOILHT7A.460
DHSL(I,N)=DHSL(I,N)-(HCAP(I)*DZ(N)+HCAPW*SMCLU0(I,N) SOILHT7A.461
& +HCAPI*SMCLF0(I,N))*DTSL(I,N) SOILHT7A.462
& -DSMCLF(I,N)*((HCAPI-HCAPW)*TSL0(I,N)-LF) SOILHT7A.463
SOILHT7A.464
ENDDO SOILHT7A.465
ENDDO SOILHT7A.466
ENDDO SOILHT7A.467
SOILHT7A.468
!----------------------------------------------------------------------- SOILHT7A.469
! Diagnose soil temperatures (K) and fractional values of unfrozen and SOILHT7A.470
! frozen water. SOILHT7A.471
!----------------------------------------------------------------------- SOILHT7A.472
SOILHT7A.473
DO N=1,NSHYD SOILHT7A.474
DO J=1,SOIL_PTS SOILHT7A.475
I=SOIL_INDEX(J) SOILHT7A.476
TSOIL(I,N)=TSL(I,N)+ZERODEGC SOILHT7A.477
STHU(I,N)=SMCLU(I,N)/SMCLSAT(I,N) SOILHT7A.478
STHF(I,N)=SMCLF(I,N)/SMCLSAT(I,N) SOILHT7A.479
ENDDO SOILHT7A.480
ENDDO SOILHT7A.481
IF (LTIMER) THEN SOILHT7A.482
CALL TIMER('SOILHTC ',104) GPB8F405.42
ENDIF SOILHT7A.484
SOILHT7A.485
RETURN SOILHT7A.486
END SOILHT7A.487
*ENDIF SOILHT7A.488