phylmd/Interface_surf/soil.f

module soil_m

  IMPLICIT NONE

contains

  SUBROUTINE soil(dtime, nisurf, snow, tsurf, tsoil, soilcap, soilflux)

    ! From LMDZ4/libf/phylmd/soil.F, version 1.1.1.1, 2004/05/19

    ! Author: Frederic Hourdin, January 30th, 1992

    ! Object: computation of the soil temperature evolution, the heat
    ! capacity per unit surface and the surface conduction flux

    ! Method: implicit time integration

    ! Consecutive ground temperatures are related by:
    ! T(k + 1) = C(k) + D(k) * T(k) (1)
    ! The coefficients C and D are computed at the t - dt time-step.
    ! Structure of the procedure:
    ! 1) new temperatures are computed using (1)
    ! 2) C and D coefficients are computed from the new temperature
    ! profile for the t + dt time-step
    ! 3) the coefficients A and B are computed where the diffusive
    ! fluxes at the t + dt time-step is given by
    ! Fdiff = A + B Ts(t + dt)
    ! or 
    ! Fdiff = F0 + Soilcap (Ts(t + dt) - Ts(t)) / dt
    ! with 
    ! F0 = A + B (Ts(t))
    ! Soilcap = B * dt

    USE indicesol, only: nbsrf, is_lic, is_oce, is_sic, is_ter
    USE dimphy, only: klon
    USE dimsoil, only: nsoilmx
    USE suphec_m, only: rtt

    REAL, intent(in):: dtime ! physical timestep (s)
    INTEGER, intent(in):: nisurf ! sub-surface index
    REAL, intent(in):: snow(:) ! (knon)
    REAL, intent(in):: tsurf(:) ! (knon) surface temperature at time-step t (K)

    real, intent(inout):: tsoil(:, :) ! (knon, nsoilmx)
    ! temperature inside the ground (K)

    REAL, intent(out):: soilcap(:) ! (knon)
    ! specific heat per unit surface (W m-2 s K-1)

    REAL, intent(out):: soilflux(:) ! (knon) 
    ! surface diffusive flux from ground (W m-2)

    ! Local:

    INTEGER knon, ig, jk
    REAL zdz2(nsoilmx)
    real z1(size(tsurf), nbsrf) ! (knon, nbsrf)
    REAL min_period, dalph_soil
    REAL ztherm_i(size(tsurf)) ! (knon)
    REAL, save:: dz1(nsoilmx), dz2(nsoilmx)
    REAL, save:: zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf)
    REAL, save:: lambda
    LOGICAL:: firstsurf(nbsrf) = .TRUE.
    REAL:: isol = 2000., isno = 2000., iice = 2000.

    ! Depths:
    REAL rk, fz1, rk1, rk2

    !-----------------------------------------------------------------------

    knon = size(tsurf)

    ! Calcul de l'inertie thermique. On initialise \`a iice m\^eme
    ! au-dessus d'un point de mer au cas o\`u le point de mer devienne
    ! point de glace au pas suivant. On corrige si on a un point de
    ! terre avec ou sans glace.

    IF (nisurf==is_sic) THEN
       DO ig = 1, knon
          ztherm_i(ig) = iice
          IF (snow(ig) > 0.0) ztherm_i(ig) = isno
       END DO
    ELSE IF (nisurf==is_lic) THEN
       DO ig = 1, knon
          ztherm_i(ig) = iice
          IF (snow(ig) > 0.0) ztherm_i(ig) = isno
       END DO
    ELSE IF (nisurf==is_ter) THEN
       DO ig = 1, knon
          ztherm_i(ig) = isol
          IF (snow(ig) > 0.0) ztherm_i(ig) = isno
       END DO
    ELSE IF (nisurf==is_oce) THEN
       DO ig = 1, knon
          ztherm_i(ig) = iice
       END DO
    ELSE
       PRINT *, 'valeur d indice non prevue', nisurf
       STOP 1
    END IF

    IF (firstsurf(nisurf)) THEN
       ! ground levels
       ! grnd=z / l where l is the skin depth of the diurnal cycle:

       min_period = 1800. ! en secondes
       dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.

       OPEN(99, FILE='soil.def', STATUS='old', FORM='formatted', ERR=9999)
       READ(99, *) min_period
       READ(99, *) dalph_soil
       PRINT *, 'Discretization for the soil model'
       PRINT *, 'First level e-folding depth', min_period, ' dalph', &
            dalph_soil
       CLOSE(99)
9999   CONTINUE

       ! la premiere couche represente un dixieme de cycle diurne
       fz1 = sqrt(min_period / 3.14)

       DO jk = 1, nsoilmx
          rk1 = jk
          rk2 = jk - 1
          dz2(jk) = fz(rk1) - fz(rk2)
       END DO
       DO jk = 1, nsoilmx - 1
          rk1 = jk + .5
          rk2 = jk - .5
          dz1(jk) = 1. / (fz(rk1) - fz(rk2))
       END DO
       lambda = fz(.5) * dz1(1)
       PRINT *, 'full layers, intermediate layers (seconds)'
       DO jk = 1, nsoilmx
          rk = jk
          rk1 = jk + .5
          rk2 = jk - .5
          PRINT *, 'fz=', fz(rk1) * fz(rk2) * 3.14, fz(rk) * fz(rk) * 3.14
       END DO
       ! PB
       firstsurf(nisurf) = .FALSE.
    ELSE
       ! Computation of the soil temperatures using the Zc and Zd
       ! coefficient computed at the previous time-step:

       ! surface temperature
       DO ig = 1, knon
          tsoil(ig, 1) = (lambda * zc(ig, 1, nisurf) + tsurf(ig)) &
               / (lambda * (1. - zd(ig, 1, nisurf)) + 1.)
       END DO

       ! other temperatures
       DO jk = 1, nsoilmx - 1
          DO ig = 1, knon
             tsoil(ig, jk + 1) = zc(ig, jk, nisurf) &
                  + zd(ig, jk, nisurf) * tsoil(ig, jk)
          END DO
       END DO
    END IF

    ! Computation of the Zc and Zd coefficient for the next step:

    IF (nisurf==is_sic) THEN
       DO ig = 1, knon
          tsoil(ig, nsoilmx) = rtt - 1.8
       END DO
    END IF

    DO jk = 1, nsoilmx
       zdz2(jk) = dz2(jk) / dtime
    END DO

    DO ig = 1, knon
       z1(ig, nisurf) = zdz2(nsoilmx) + dz1(nsoilmx - 1)
       zc(ig, nsoilmx - 1, nisurf) = zdz2(nsoilmx) * tsoil(ig, nsoilmx) / &
            z1(ig, nisurf)
       zd(ig, nsoilmx - 1, nisurf) = dz1(nsoilmx - 1) / z1(ig, nisurf)
    END DO

    DO jk = nsoilmx - 1, 2, - 1
       DO ig = 1, knon
          z1(ig, nisurf) = 1. / (zdz2(jk) + dz1(jk - 1) &
               + dz1(jk) * (1. - zd(ig, jk, nisurf)))
          zc(ig, jk - 1, nisurf) = (tsoil(ig, jk) * zdz2(jk) &
               + dz1(jk) * zc(ig, jk, nisurf)) * z1(ig, nisurf)
          zd(ig, jk - 1, nisurf) = dz1(jk - 1) * z1(ig, nisurf)
       END DO
    END DO

    ! computation of the surface diffusive flux from ground and
    ! calorific capacity of the ground:

    DO ig = 1, knon
       soilflux(ig) = ztherm_i(ig) * dz1(1) * (zc(ig, 1, nisurf) + (zd(ig, 1, &
            nisurf) - 1.) * tsoil(ig, 1))
       soilcap(ig) = ztherm_i(ig) * (dz2(1) &
            + dtime * (1. - zd(ig, 1, nisurf)) * dz1(1))
       z1(ig, nisurf) = lambda * (1. - zd(ig, 1, nisurf)) + 1.
       soilcap(ig) = soilcap(ig) / z1(ig, nisurf)
       soilflux(ig) = soilflux(ig) + soilcap(ig) * (tsoil(ig, 1) &
            * z1(ig, nisurf) - lambda * zc(ig, 1, nisurf) - tsurf(ig)) / dtime
    END DO

  contains

    pure real function fz(rk)

      real, intent(in):: rk

      !-----------------------------------------

      fz = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)

    end function fz

  END SUBROUTINE soil

end module soil_m
1	guez	101	module soil_m
2	guez	3
3	guez	81	IMPLICIT NONE
4	guez	3
5	guez	101	contains
6	guez	3
7	guez	207	SUBROUTINE soil(dtime, nisurf, snow, tsurf, tsoil, soilcap, soilflux)
8	guez	3
9	guez	207	! From LMDZ4/libf/phylmd/soil.F, version 1.1.1.1, 2004/05/19
10	guez	3
11	guez	207	! Author: Frederic Hourdin, January 30th, 1992
12	guez	3
13	guez	207	! Object: computation of the soil temperature evolution, the heat
14			! capacity per unit surface and the surface conduction flux
15	guez	3
16	guez	175	! Method: implicit time integration
17	guez	3
18	guez	101	! Consecutive ground temperatures are related by:
19	guez	207	! T(k + 1) = C(k) + D(k) * T(k) (1)
20			! The coefficients C and D are computed at the t - dt time-step.
21			! Structure of the procedure:
22	guez	175	! 1) new temperatures are computed using (1)
23			! 2) C and D coefficients are computed from the new temperature
24	guez	207	! profile for the t + dt time-step
25	guez	175	! 3) the coefficients A and B are computed where the diffusive
26	guez	207	! fluxes at the t + dt time-step is given by
27			! Fdiff = A + B Ts(t + dt)
28			! or
29			! Fdiff = F0 + Soilcap (Ts(t + dt) - Ts(t)) / dt
30			! with
31			! F0 = A + B (Ts(t))
32			! Soilcap = B * dt
33	guez	3
34	guez	207	USE indicesol, only: nbsrf, is_lic, is_oce, is_sic, is_ter
35			USE dimphy, only: klon
36			USE dimsoil, only: nsoilmx
37			USE suphec_m, only: rtt
38	guez	175
39	guez	207	REAL, intent(in):: dtime ! physical timestep (s)
40			INTEGER, intent(in):: nisurf ! sub-surface index
41			REAL, intent(in):: snow(:) ! (knon)
42			REAL, intent(in):: tsurf(:) ! (knon) surface temperature at time-step t (K)
43	guez	3
44	guez	207	real, intent(inout):: tsoil(:, :) ! (knon, nsoilmx)
45			! temperature inside the ground (K)
46	guez	3
47	guez	207	REAL, intent(out):: soilcap(:) ! (knon)
48			! specific heat per unit surface (W m-2 s K-1)
49	guez	3
50	guez	207	REAL, intent(out):: soilflux(:) ! (knon)
51			! surface diffusive flux from ground (W m-2)
52	guez	3
53	guez	207	! Local:
54	guez	3
55	guez	207	INTEGER knon, ig, jk
56			REAL zdz2(nsoilmx)
57			real z1(size(tsurf), nbsrf) ! (knon, nbsrf)
58	guez	101	REAL min_period, dalph_soil
59	guez	207	REAL ztherm_i(size(tsurf)) ! (knon)
60			REAL, save:: dz1(nsoilmx), dz2(nsoilmx)
61			REAL, save:: zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf)
62			REAL, save:: lambda
63			LOGICAL:: firstsurf(nbsrf) = .TRUE.
64			REAL:: isol = 2000., isno = 2000., iice = 2000.
65	guez	3
66	guez	207	! Depths:
67			REAL rk, fz1, rk1, rk2
68	guez	3
69	guez	207	!-----------------------------------------------------------------------
70	guez	3
71	guez	207	knon = size(tsurf)
72	guez	3
73	guez	207	! Calcul de l'inertie thermique. On initialise \`a iice m\^eme
74			! au-dessus d'un point de mer au cas o\`u le point de mer devienne
75			! point de glace au pas suivant. On corrige si on a un point de
76			! terre avec ou sans glace.
77	guez	3
78	guez	175	IF (nisurf==is_sic) THEN
79	guez	101	DO ig = 1, knon
80			ztherm_i(ig) = iice
81	guez	207	IF (snow(ig) > 0.0) ztherm_i(ig) = isno
82	guez	101	END DO
83	guez	175	ELSE IF (nisurf==is_lic) THEN
84	guez	101	DO ig = 1, knon
85			ztherm_i(ig) = iice
86	guez	207	IF (snow(ig) > 0.0) ztherm_i(ig) = isno
87	guez	101	END DO
88	guez	175	ELSE IF (nisurf==is_ter) THEN
89	guez	101	DO ig = 1, knon
90			ztherm_i(ig) = isol
91	guez	207	IF (snow(ig) > 0.0) ztherm_i(ig) = isno
92	guez	101	END DO
93	guez	175	ELSE IF (nisurf==is_oce) THEN
94	guez	101	DO ig = 1, knon
95			ztherm_i(ig) = iice
96			END DO
97			ELSE
98	guez	175	PRINT *, 'valeur d indice non prevue', nisurf
99	guez	101	STOP 1
100			END IF
101	guez	3
102	guez	175	IF (firstsurf(nisurf)) THEN
103	guez	101	! ground levels
104	guez	207	! grnd=z / l where l is the skin depth of the diurnal cycle:
105	guez	3
106	guez	101	min_period = 1800. ! en secondes
107			dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
108	guez	81
109	guez	207	OPEN(99, FILE='soil.def', STATUS='old', FORM='formatted', ERR=9999)
110			READ(99, *) min_period
111			READ(99, *) dalph_soil
112	guez	101	PRINT *, 'Discretization for the soil model'
113	guez	175	PRINT *, 'First level e-folding depth', min_period, ' dalph', &
114	guez	101	dalph_soil
115	guez	207	CLOSE(99)
116	guez	101	9999 CONTINUE
117
118			! la premiere couche represente un dixieme de cycle diurne
119	guez	207	fz1 = sqrt(min_period / 3.14)
120	guez	101
121			DO jk = 1, nsoilmx
122			rk1 = jk
123			rk2 = jk - 1
124			dz2(jk) = fz(rk1) - fz(rk2)
125			END DO
126			DO jk = 1, nsoilmx - 1
127			rk1 = jk + .5
128			rk2 = jk - .5
129	guez	207	dz1(jk) = 1. / (fz(rk1) - fz(rk2))
130	guez	101	END DO
131	guez	207	lambda = fz(.5) * dz1(1)
132	guez	101	PRINT *, 'full layers, intermediate layers (seconds)'
133			DO jk = 1, nsoilmx
134			rk = jk
135			rk1 = jk + .5
136			rk2 = jk - .5
137	guez	207	PRINT , 'fz=', fz(rk1) fz(rk2) * 3.14, fz(rk) * fz(rk) * 3.14
138	guez	101	END DO
139			! PB
140	guez	175	firstsurf(nisurf) = .FALSE.
141	guez	157	ELSE
142	guez	221	! Computation of the soil temperatures using the Zc and Zd
143	guez	101	! coefficient computed at the previous time-step:
144
145			! surface temperature
146			DO ig = 1, knon
147	guez	207	tsoil(ig, 1) = (lambda * zc(ig, 1, nisurf) + tsurf(ig)) &
148			/ (lambda * (1. - zd(ig, 1, nisurf)) + 1.)
149	guez	101	END DO
150
151			! other temperatures
152			DO jk = 1, nsoilmx - 1
153			DO ig = 1, knon
154	guez	207	tsoil(ig, jk + 1) = zc(ig, jk, nisurf) &
155			+ zd(ig, jk, nisurf) * tsoil(ig, jk)
156	guez	101	END DO
157			END DO
158	guez	207	END IF
159	guez	101
160	guez	221	! Computation of the Zc and Zd coefficient for the next step:
161	guez	81
162	guez	175	IF (nisurf==is_sic) THEN
163	guez	101	DO ig = 1, knon
164	guez	175	tsoil(ig, nsoilmx) = rtt - 1.8
165	guez	101	END DO
166			END IF
167	guez	81
168	guez	101	DO jk = 1, nsoilmx
169	guez	207	zdz2(jk) = dz2(jk) / dtime
170	guez	81	END DO
171
172			DO ig = 1, knon
173	guez	207	z1(ig, nisurf) = zdz2(nsoilmx) + dz1(nsoilmx - 1)
174			zc(ig, nsoilmx - 1, nisurf) = zdz2(nsoilmx) * tsoil(ig, nsoilmx) / &
175	guez	175	z1(ig, nisurf)
176	guez	207	zd(ig, nsoilmx - 1, nisurf) = dz1(nsoilmx - 1) / z1(ig, nisurf)
177	guez	81	END DO
178
179	guez	207	DO jk = nsoilmx - 1, 2, - 1
180	guez	101	DO ig = 1, knon
181	guez	207	z1(ig, nisurf) = 1. / (zdz2(jk) + dz1(jk - 1) &
182			+ dz1(jk) * (1. - zd(ig, jk, nisurf)))
183			zc(ig, jk - 1, nisurf) = (tsoil(ig, jk) * zdz2(jk) &
184			+ dz1(jk) * zc(ig, jk, nisurf)) * z1(ig, nisurf)
185			zd(ig, jk - 1, nisurf) = dz1(jk - 1) * z1(ig, nisurf)
186	guez	101	END DO
187			END DO
188	guez	81
189	guez	101	! computation of the surface diffusive flux from ground and
190			! calorific capacity of the ground:
191	guez	81
192			DO ig = 1, knon
193	guez	207	soilflux(ig) = ztherm_i(ig) * dz1(1) * (zc(ig, 1, nisurf) + (zd(ig, 1, &
194			nisurf) - 1.) * tsoil(ig, 1))
195			soilcap(ig) = ztherm_i(ig) * (dz2(1) &
196			+ dtime * (1. - zd(ig, 1, nisurf)) * dz1(1))
197			z1(ig, nisurf) = lambda * (1. - zd(ig, 1, nisurf)) + 1.
198			soilcap(ig) = soilcap(ig) / z1(ig, nisurf)
199			soilflux(ig) = soilflux(ig) + soilcap(ig) * (tsoil(ig, 1) &
200			* z1(ig, nisurf) - lambda * zc(ig, 1, nisurf) - tsurf(ig)) / dtime
201	guez	81	END DO
202
203	guez	108	contains
204
205	guez	207	pure real function fz(rk)
206	guez	202
207			real, intent(in):: rk
208
209			!-----------------------------------------
210
211	guez	207	fz = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
212	guez	202
213	guez	108	end function fz
214
215	guez	101	END SUBROUTINE soil
216	guez	81
217	guez	101	end module soil_m