# if defined O_green_australia || O_green_africa
subroutine MTLM_STATE (LAND_PTS, LAND_INDEX, DZ_SOIL
# else
subroutine MTLM_STATE (POINTS, LAND_PTS, LAND_INDEX, DZ_SOIL
# endif
&, HCAP_SOIL, KS, THETA_SAT, LF, TM, TIMESTEP
&, G, RAIN, SNOW, E, ESUB, M, LYING_SNOW, TS1
&, RUNOFF, SNOWMELT, MNEG)
#if defined O_mtlm
!-----------------------------------------------------------------------
! Routine to update land surface prognostic variables
! (soil moisture, soil temperature, lying snow mass).
! Also diagnoses runoff and snowmelt
!**********************************************************************
! this file is based on code that may have had the following copyright:
! (c) CROWN COPYRIGHT 1997, U.K. METEOROLOGICAL OFFICE.
! Permission has been granted by the authors to the public to copy
! and use this software without charge, provided that this Notice and
! any statement of authorship are reproduced on all copies. Neither the
! Crown nor the U.K. Meteorological Office makes any warranty, express
! or implied, or assumes any liability or responsibility for the use of
! this software.
!**********************************************************************
!-----------------------------------------------------------------------
implicit none
# if defined O_green_australia || O_green_africa
include "size.h"
include "csbc.h"
# endif
! POINTS = IN Maximum number of land points.
! LAND_PTS = IN Number of land points.
! LAND_INDEX = IN Indices of landpoints
# if defined O_green_australia || O_green_africa
integer LAND_PTS, LAND_INDEX(POINTS)
integer I, J, L
real GREEN_THETA_SAT
# else
integer POINTS, LAND_PTS, LAND_INDEX(POINTS)
integer I, L
# endif
! Surface parameters
! DZ_SOIL = IN Soil layer thickness (m).
! HCAP_SOIL = IN Soil heat capacity (W/m3/K).
! TIMESTEP = IN Timestep (s).
! KS = IN Saturated hydraulic conductivity
! THETA_SAT = IN Saturated volumetric soil moisture concentration
! (m3/m3).
! LF = Latent heat of fusion (J/kg).
! TM = Melting point of fresh water (K).
! G = IN Ground heat flux (W/m2).
real DZ_SOIL, HCAP_SOIL, TIMESTEP, KS, THETA_SAT, LF, TM
real G(POINTS)
! Driving variables
! RAIN = IN Rainfall (kg/m2/s).
! SNOW = IN Snowfall (kg/m2/s).
! E = IN Evapotranspiration (kg/m2/s).
! ESUB = IN Sublimation (kg/m2/s).
! M = INOUT Soil moisture (kg/m2).
! NEG = INOUT Negative Soil moisture (kg/m2).
! LYING_SNOW = INOUT Lying snow (kg/m2).
! TS1 = INOUT Soil temperature (K).
! RUNOFF = OUT Runoff (kg/m2/s).
! SNOWMELT = OUT Snow melt (kg/m2/s).
real RAIN(POINTS), SNOW(POINTS), E(POINTS), ESUB(POINTS)
real M(POINTS), MNEG(POINTS), LYING_SNOW(POINTS)
real TS1(POINTS), RUNOFF(POINTS), SNOWMELT(POINTS)
! Local parameters
! B = Clapp-Hornberger exponent.
real B
parameter (B=6.6)
!----------------------------------------------------------------------
! Update the soil temperature and diagnose snowmelt.
!----------------------------------------------------------------------
!CDIR NODEP
do I=1,LAND_PTS
L = LAND_INDEX(I)
TS1(L) = TS1(L) + TIMESTEP* G(L) / (DZ_SOIL*HCAP_SOIL)
SNOWMELT(L) = 0.0
if ((LYING_SNOW(L).gt.0.0) .and. (TS1(L).gt.TM)) then
SNOWMELT(L) = HCAP_SOIL*DZ_SOIL*(TS1(L)-TM)
& / (LF*TIMESTEP)
if ((SNOWMELT(L)).gt.(SNOW(L)-ESUB(L)+LYING_SNOW(L)/TIMESTEP))
& then
! limit snowmelt to amount of snow and fix soil temperature
SNOWMELT(L) = SNOW(L)-ESUB(L)+LYING_SNOW(L)/TIMESTEP
TS1(L) = TS1(L) - SNOWMELT(L)*LF*TIMESTEP/HCAP_SOIL*DZ_SOIL
else
TS1(L) = TM
endif
endif
!----------------------------------------------------------------------
! Update the lying snow.
!----------------------------------------------------------------------
LYING_SNOW(L) = LYING_SNOW(L)
& + TIMESTEP*(SNOW(L)-ESUB(L)-SNOWMELT(L))
if (LYING_SNOW(L) .lt. 0.) then
! if negative snow, change extra sublimation to evaporation
ESUB(L) = ESUB(L) + LYING_SNOW(L)/TIMESTEP
E(L) = E(L) - LYING_SNOW(L)/TIMESTEP
TS1(L) = TS1(L) + LF*LYING_SNOW(L)/(DZ_SOIL*HCAP_SOIL)
LYING_SNOW(L) = 0.
endif
!----------------------------------------------------------------------
! Update the soil moisture.
!----------------------------------------------------------------------
RUNOFF(L) = KS*(M(L)/(1000.0*DZ_SOIL*THETA_SAT))**(2*B+3)
M(L) = M(L) + TIMESTEP*(RAIN(L)+SNOWMELT(L)-E(L)-RUNOFF(L))
! keep track of any negative moisture for conservation
if (M(L) + MNEG(L) .lt. 0.) then
MNEG(L) = MNEG(L) + M(L)
M(L) = 0.
else
M(L) = MNEG(L) + M(L)
MNEG(L) = 0.
endif
enddo
# if defined O_green_australia
GREEN_THETA_SAT = 0.35
! Irrigate Australia and watch the plants grow
do J=29,44
do I=33,43
L = land_map(I,J)
M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
enddo
enddo
# endif
# if defined O_green_africa
GREEN_THETA_SAT = 0.35
! Irrigate N. Africa and watch the plants grow
! VSAT = 0.458m3 H2O/m3 soil = theta_sat
! RHO_WATER = 1000kg/m3
! ROOTDEP = dz_soil
! First do the Western part
do J=61,72
do I=96,101
L = land_map(I,J)
M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
enddo
enddo
! Then most of the Eastern part
do J=61,72
do I=1,9
L = land_map(I,J)
M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
enddo
enddo
! Now fill in an area to the south
J=60
do I=4,11
L = land_map(I,J)
M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
enddo
! Now fill in more Eastern area (Eritrea and Sudan)
do J=61,69
do I=10,11
L = land_map(I,J)
M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
enddo
enddo
! Now fill in more Eastern area (Egypt and Sudan)
do J=60,64
I=12
L = land_map(I,J)
M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
enddo
# if defined O_green_Saudi_Arabia
!
do J=60,68
do I=12,14
L = land_map(I,J)
M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
enddo
enddo
!
do J=60,68
I=15
L = land_map(I,J)
M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
enddo
do J=61,66
do I=16,17
L = land_map(I,J)
M(L) = GREEN_THETA_SAT*1000.*DZ_SOIL
enddo
enddo
# endif
# endif
#endif
return
end