SUBROUTINE SOIL_HYD ( 3,4SOILHY5A.23
& NPNTS,NSHYD,SOIL_PTS,SOIL_INDEX,B,DZ,EXT,FW,KS,SATHH,TIMESTEP SOILHY5A.24
&,V_SAT,SLOW_RUNOFF,SMCL,STHU,W_FLUX,STF_SLOW_RUNOFF,LTIMER SOILHY5A.25
&) SOILHY5A.26
SOILHY5A.27
IMPLICIT NONE SOILHY5A.28
! SOILHY5A.29
! Description: SOILHY5A.30
! Increments the layer soil moisture contents and calculates SOILHY5A.31
! calculates gravitational runoff. Calls the following: SOILHY5A.32
! SOILHY5A.33
! HYD_CON - to calculate the hydraulic conductivity SOILHY5A.34
! (Cox, 6/95) SOILHY5A.35
! SOILHY5A.36
! DARCY - to calculate the Darcian fluxes between soil layers SOILHY5A.37
! (Cox, 6/95) SOILHY5A.38
! SOILHY5A.39
! SOILHY5A.40
! Documentation : UM Documentation Paper 25 SOILHY5A.41
! SOILHY5A.42
! Current Code Owner : David Gregory SOILHY5A.43
! SOILHY5A.44
! History: SOILHY5A.45
! Version Date Comment SOILHY5A.46
! ------- ---- ------- SOILHY5A.47
! 4.1 6/96 New deck. Peter Cox SOILHY5A.48
! SOILHY5A.49
! Code Description: SOILHY5A.50
! Language: FORTRAN 77 + common extensions. SOILHY5A.51
! SOILHY5A.52
! System component covered: P25 SOILHY5A.53
! System Task: P25 SOILHY5A.54
! SOILHY5A.55
SOILHY5A.56
! Global variables: SOILHY5A.57
*CALL C_DENSTY 
SOILHY5A.58
*CALL C_POND SOILHY5A.59
*CALL C_SOILH SOILHY5A.60
SOILHY5A.61
! Subroutine arguments SOILHY5A.62
! Scalar arguments with intent(IN) : SOILHY5A.63
INTEGER SOILHY5A.64
& NPNTS ! IN Number of gridpoints. SOILHY5A.65
&,NSHYD ! IN Number of soil moisture levels. SOILHY5A.66
&,SOIL_PTS ! IN Number of soil points. SOILHY5A.67
SOILHY5A.68
REAL SOILHY5A.69
& TIMESTEP ! IN Model timestep (s). SOILHY5A.70
SOILHY5A.71
LOGICAL SOILHY5A.72
& STF_SLOW_RUNOFF ! IN Stash flag for sub-surface runoff. SOILHY5A.73
SOILHY5A.74
SOILHY5A.75
SOILHY5A.76
! Array arguments with intent(IN) : SOILHY5A.77
INTEGER SOILHY5A.78
& SOIL_INDEX(NPNTS) ! IN Array of soil points. SOILHY5A.79
SOILHY5A.80
REAL SOILHY5A.81
& B(NPNTS) ! IN Clapp-Hornberger exponent. SOILHY5A.82
&,DZ(NSHYD) ! IN Thicknesses of the soil layers (m). SOILHY5A.83
&,EXT(NPNTS,NSHYD) ! IN Extraction of water from each soil SOILHY5A.84
! ! layer (kg/m2/s). SOILHY5A.85
&,FW(NPNTS) ! IN Throughfall from canopy plus snowmelt SOILHY5A.86
! ! minus surface runoff (kg/m2/s). SOILHY5A.87
&,KS(NPNTS) ! IN Saturated hydraulic conductivity SOILHY5A.88
! ! (kg/m2/s). SOILHY5A.89
&,SATHH(NPNTS) ! IN Saturated soil water pressure (m). SOILHY5A.90
&,V_SAT(NPNTS) ! IN Volumetric soil moisture SOILHY5A.91
! ! concentration at saturation SOILHY5A.92
! ! (m3 H2O/m3 soil). SOILHY5A.93
! SOILHY5A.94
LOGICAL LTIMER ! Logical switch for TIMER diags SOILHY5A.95
SOILHY5A.96
! Array arguments with intent(OUT) : SOILHY5A.97
REAL SOILHY5A.98
& SLOW_RUNOFF(NPNTS) ! OUT Drainage from the base of the SOILHY5A.99
! ! soil profile (kg/m2/s). SOILHY5A.100
&,SMCLSAT(NPNTS,NSHYD) ! OUT The saturation moisture content of SOILHY5A.101
! ! each layer (kg/m2). SOILHY5A.102
&,W_FLUX(NPNTS,0:NSHYD)! OUT The fluxes of water between layers SOILHY5A.103
! ! (kg/m2/s). SOILHY5A.104
SOILHY5A.105
! Array arguments with intent(INOUT) : SOILHY5A.106
REAL SOILHY5A.107
& SMCL(NPNTS,NSHYD) ! INOUT Total soil moisture contents SOILHY5A.108
! ! of each layer (kg/m2). SOILHY5A.109
&,STHU(NPNTS,NSHYD) ! INOUT Unfrozen soil moisture content of ea SOILHY5A.110
! ! layer as a fraction of saturation. SOILHY5A.111
SOILHY5A.112
SOILHY5A.113
! Local scalars: SOILHY5A.114
INTEGER SOILHY5A.115
& I,J,N ! WORK Loop counters. SOILHY5A.116
SOILHY5A.117
! Local arrays: SOILHY5A.118
REAL SOILHY5A.119
& DSMCL(NPNTS) ! WORK The transfer of soil SOILHY5A.120
! ! moisture (kg/m2/timestep). SOILHY5A.121
&,EXCESS(NPNTS) ! WORK Excess soil moisture (kg/m2). SOILHY5A.122
&,SMCLMAX(NPNTS,NSHYD) ! WORK The maximum moisture content SOILHY5A.123
! ! of each layer (kg/m2). SOILHY5A.124
&,SMCLU(NPNTS,NSHYD) ! WORK Unfrozen soil moisture contents SOILHY5A.125
! ! of each layer (kg/m2). SOILHY5A.126
&,STHUK(NPNTS) ! WORK Fractional saturation of lowest SOILHY5A.127
! ! layer. SOILHY5A.128
SOILHY5A.129
! Function & Subroutine calls: SOILHY5A.130
EXTERNAL SOILHY5A.131
& HYD_CON,DARCY SOILHY5A.132
SOILHY5A.133
IF (LTIMER) THEN SOILHY5A.134
CALL TIMER('SOILHYD ',103) GPB8F405.150
ENDIF SOILHY5A.136
!----------------------------------------------------------------------- SOILHY5A.137
! Calculate the unfrozen soil moisture contents and the saturation SOILHY5A.138
! total soil moisture for each layer. SOILHY5A.139
!----------------------------------------------------------------------- SOILHY5A.140
! CDIR$ IVDEP here would force vectorization but changes results! SOILHY5A.141
! Fujitsu vectorization directive GRB0F405.537
!OCL NOVREC GRB0F405.538
DO N=1,NSHYD SOILHY5A.142
DO J=1,SOIL_PTS SOILHY5A.143
I=SOIL_INDEX(J) SOILHY5A.144
SMCLSAT(I,N)=RHO_WATER*DZ(N)*V_SAT(I) SOILHY5A.145
SMCLU(I,N)=STHU(I,N)*SMCLSAT(I,N) SOILHY5A.146
ENDDO SOILHY5A.147
ENDDO SOILHY5A.148
SOILHY5A.149
!----------------------------------------------------------------------- SOILHY5A.150
! Top boundary condition SOILHY5A.151
!----------------------------------------------------------------------- SOILHY5A.152
DO J=1,SOIL_PTS SOILHY5A.153
I=SOIL_INDEX(J) SOILHY5A.154
W_FLUX(I,0)=FW(I) SOILHY5A.155
ENDDO SOILHY5A.156
SOILHY5A.157
!----------------------------------------------------------------------- SOILHY5A.158
! Define the maximum water content that may exist in each layer SOILHY5A.159
!----------------------------------------------------------------------- SOILHY5A.160
DO N=2,NSHYD SOILHY5A.161
DO J=1,SOIL_PTS SOILHY5A.162
I=SOIL_INDEX(J) SOILHY5A.163
SMCLMAX(I,N)=SMCLSAT(I,N) SOILHY5A.164
ENDDO SOILHY5A.165
ENDDO SOILHY5A.166
SOILHY5A.167
!----------------------------------------------------------------------- SOILHY5A.168
! Allow for some ponding of water by permitting excess moisture to SOILHY5A.169
! remain in the top layer. SOILHY5A.170
!----------------------------------------------------------------------- SOILHY5A.171
DO J=1,SOIL_PTS SOILHY5A.172
I=SOIL_INDEX(J) SOILHY5A.173
SMCLMAX(I,1)=SMCLSAT(I,1)+POND_MAX SOILHY5A.174
ENDDO SOILHY5A.175
SOILHY5A.176
!----------------------------------------------------------------------- SOILHY5A.177
! Increment the soil moisture contents for each layer, beginning with SOILHY5A.178
! the bottom. SOILHY5A.179
!----------------------------------------------------------------------- SOILHY5A.180
DO J=1,SOIL_PTS SOILHY5A.181
I=SOIL_INDEX(J) SOILHY5A.182
STHUK(I)=STHU(I,NSHYD) SOILHY5A.183
ENDDO SOILHY5A.184
CALL HYD_CON (NPNTS,SOIL_PTS,SOIL_INDEX,B,KS,STHUK, SOILHY5A.185
& W_FLUX(1,NSHYD),LTIMER) SOILHY5A.186
SOILHY5A.187
CDIR$ IVDEP SOILHY5A.188
! Fujitsu vectorization directive GRB0F405.539
!OCL NOVREC GRB0F405.540
DO J=1,SOIL_PTS SOILHY5A.189
I=SOIL_INDEX(J) SOILHY5A.190
DSMCL(I)=(W_FLUX(I,NSHYD)+EXT(I,NSHYD))*TIMESTEP SOILHY5A.191
IF (DSMCL(I).GT.SMCLU(I,NSHYD)) THEN SOILHY5A.192
DSMCL(I)=SMCLU(I,NSHYD) SOILHY5A.193
W_FLUX(I,NSHYD)=DSMCL(I)/TIMESTEP-EXT(I,NSHYD) SOILHY5A.194
ENDIF SOILHY5A.195
SOILHY5A.196
SMCL(I,NSHYD)=SMCL(I,NSHYD)-TIMESTEP*(W_FLUX(I,NSHYD) SOILHY5A.197
& +EXT(I,NSHYD)) SOILHY5A.198
SMCLU(I,NSHYD)=SMCLU(I,NSHYD)-TIMESTEP*(W_FLUX(I,NSHYD) SOILHY5A.199
& +EXT(I,NSHYD)) SOILHY5A.200
STHU(I,NSHYD)=SMCLU(I,NSHYD)/SMCLSAT(I,NSHYD) SOILHY5A.201
ENDDO SOILHY5A.202
SOILHY5A.203
!----------------------------------------------------------------------- SOILHY5A.204
! Layers NSHYD to 1 SOILHY5A.205
!----------------------------------------------------------------------- SOILHY5A.206
DO N=NSHYD,2,-1 SOILHY5A.207
CALL DARCY (NPNTS,SOIL_PTS,SOIL_INDEX,B,KS,SATHH SOILHY5A.208
&, STHU(1,N-1),DZ(N-1),STHU(1,N),DZ(N) SOILHY5A.209
&, W_FLUX(1,N-1),LTIMER) SOILHY5A.210
SOILHY5A.211
CDIR$ IVDEP SOILHY5A.212
! Fujitsu vectorization directive GRB0F405.541
!OCL NOVREC GRB0F405.542
DO J=1,SOIL_PTS SOILHY5A.213
I=SOIL_INDEX(J) SOILHY5A.214
DSMCL(I)=(W_FLUX(I,N-1)+EXT(I,N-1))*TIMESTEP SOILHY5A.215
IF (DSMCL(I).GT.SMCLU(I,N-1)) THEN SOILHY5A.216
DSMCL(I)=SMCLU(I,N-1) SOILHY5A.217
W_FLUX(I,N-1)=DSMCL(I)/TIMESTEP-EXT(I,N-1) SOILHY5A.218
ELSEIF (DSMCL(I).LT.(SMCL(I,N-1)-SMCLMAX(I,N-1))) THEN SOILHY5A.219
DSMCL(I)=SMCL(I,N-1)-SMCLMAX(I,N-1) SOILHY5A.220
W_FLUX(I,N-1)=DSMCL(I)/TIMESTEP-EXT(I,N-1) SOILHY5A.221
ENDIF SOILHY5A.222
SOILHY5A.223
DSMCL(I)=W_FLUX(I,N-1)*TIMESTEP SOILHY5A.224
IF (DSMCL(I).GT.(SMCLMAX(I,N)-SMCL(I,N))) THEN SOILHY5A.225
DSMCL(I)=SMCLMAX(I,N)-SMCL(I,N) SOILHY5A.226
W_FLUX(I,N-1)=DSMCL(I)/TIMESTEP SOILHY5A.227
ELSEIF (DSMCL(I).LT.(-SMCLU(I,N))) THEN SOILHY5A.228
DSMCL(I)=-SMCLU(I,N) SOILHY5A.229
W_FLUX(I,N-1)=DSMCL(I)/TIMESTEP SOILHY5A.230
ENDIF SOILHY5A.231
SOILHY5A.232
SMCL(I,N)=SMCL(I,N)+TIMESTEP*W_FLUX(I,N-1) SOILHY5A.233
SMCLU(I,N)=SMCLU(I,N)+TIMESTEP*W_FLUX(I,N-1) SOILHY5A.234
STHU(I,N)=SMCLU(I,N)/SMCLSAT(I,N) SOILHY5A.235
SMCL(I,N-1)=SMCL(I,N-1)-TIMESTEP*(W_FLUX(I,N-1) SOILHY5A.236
& +EXT(I,N-1)) SOILHY5A.237
SMCLU(I,N-1)=SMCLU(I,N-1)-TIMESTEP*(W_FLUX(I,N-1) SOILHY5A.238
& +EXT(I,N-1)) SOILHY5A.239
STHU(I,N-1)=SMCLU(I,N-1)/SMCLSAT(I,N-1) SOILHY5A.240
SOILHY5A.241
ENDDO SOILHY5A.242
ENDDO SOILHY5A.243
SOILHY5A.244
CDIR$ IVDEP SOILHY5A.245
! Fujitsu vectorization directive GRB0F405.543
!OCL NOVREC GRB0F405.544
DO J=1,SOIL_PTS SOILHY5A.246
I=SOIL_INDEX(J) SOILHY5A.247
SMCL(I,1)=SMCL(I,1)+TIMESTEP*W_FLUX(I,0) SOILHY5A.248
SMCLU(I,1)=SMCLU(I,1)+TIMESTEP*W_FLUX(I,0) SOILHY5A.249
STHU(I,1)=SMCLU(I,1)/SMCLSAT(I,1) SOILHY5A.250
ENDDO SOILHY5A.251
SOILHY5A.252
!----------------------------------------------------------------------- SOILHY5A.253
! If a layer is supersaturated move the excess moisture to the layer SOILHY5A.254
! below. SOILHY5A.255
!----------------------------------------------------------------------- SOILHY5A.256
DO N=1,NSHYD-1 SOILHY5A.257
CDIR$ IVDEP SOILHY5A.258
! Fujitsu vectorization directive GRB0F405.545
!OCL NOVREC GRB0F405.546
DO J=1,SOIL_PTS SOILHY5A.259
I=SOIL_INDEX(J) SOILHY5A.260
SOILHY5A.261
EXCESS(I)=MAX((SMCL(I,N)-SMCLMAX(I,N)),0.0) SOILHY5A.262
SOILHY5A.263
IF (EXCESS(I).GT.0.0) THEN SOILHY5A.264
SOILHY5A.265
W_FLUX(I,N)=W_FLUX(I,N)+EXCESS(I)/TIMESTEP SOILHY5A.266
SOILHY5A.267
SMCL(I,N)=SMCL(I,N)-EXCESS(I) SOILHY5A.268
SMCLU(I,N)=SMCLU(I,N)-EXCESS(I) SOILHY5A.269
STHU(I,N)=SMCLU(I,N)/SMCLSAT(I,N) SOILHY5A.270
SOILHY5A.271
SMCL(I,N+1)=SMCL(I,N+1)+EXCESS(I) SOILHY5A.272
SMCLU(I,N+1)=SMCLU(I,N+1)+EXCESS(I) SOILHY5A.273
STHU(I,N+1)=SMCLU(I,N+1)/SMCLSAT(I,N+1) SOILHY5A.274
SOILHY5A.275
ENDIF SOILHY5A.276
SOILHY5A.277
ENDDO SOILHY5A.278
ENDDO SOILHY5A.279
SOILHY5A.280
!----------------------------------------------------------------------- SOILHY5A.281
! If there is still excess moisture add this to the flux from the base SOILHY5A.282
! of the soil profile. SOILHY5A.283
!----------------------------------------------------------------------- SOILHY5A.284
CDIR$ IVDEP SOILHY5A.285
! Fujitsu vectorization directive GRB0F405.547
!OCL NOVREC GRB0F405.548
DO J=1,SOIL_PTS SOILHY5A.286
I=SOIL_INDEX(J) SOILHY5A.287
EXCESS(I)=MAX((SMCL(I,NSHYD)-SMCLMAX(I,NSHYD)),0.0) SOILHY5A.288
IF (EXCESS(I).GT.0.0) THEN SOILHY5A.289
SOILHY5A.290
W_FLUX(I,NSHYD)=W_FLUX(I,NSHYD)+EXCESS(I)/TIMESTEP SOILHY5A.291
SOILHY5A.292
SMCL(I,NSHYD)=SMCL(I,NSHYD)-EXCESS(I) SOILHY5A.293
SMCLU(I,NSHYD)=SMCLU(I,NSHYD)-EXCESS(I) SOILHY5A.294
STHU(I,NSHYD)=SMCLU(I,NSHYD)/SMCLSAT(I,NSHYD) SOILHY5A.295
SOILHY5A.296
ENDIF SOILHY5A.297
ENDDO SOILHY5A.298
SOILHY5A.299
!----------------------------------------------------------------------- SOILHY5A.300
! Output slow runoff (drainage) diagnostic. SOILHY5A.301
!----------------------------------------------------------------------- SOILHY5A.302
IF (STF_SLOW_RUNOFF) THEN SOILHY5A.303
DO I=1,NPNTS SOILHY5A.304
SLOW_RUNOFF(I)=0.0 SOILHY5A.305
ENDDO SOILHY5A.306
DO J=1,SOIL_PTS SOILHY5A.307
I=SOIL_INDEX(J) SOILHY5A.308
SLOW_RUNOFF(I)=W_FLUX(I,NSHYD) SOILHY5A.309
ENDDO SOILHY5A.310
ENDIF SOILHY5A.311
IF (LTIMER) THEN SOILHY5A.312
CALL TIMER('SOILHYD ',104) GPB8F405.151
ENDIF SOILHY5A.314
RETURN SOILHY5A.315
END SOILHY5A.316
*ENDIF SOILHY5A.317