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	module soil_m
2
3	IMPLICIT NONE
4
5	contains
6
7	SUBROUTINE soil(dtime, nisurf, snow, tsurf, tsoil, soilcap, soilflux)
8
9	! From LMDZ4/libf/phylmd/soil.F, version 1.1.1.1, 2004/05/19
10
11	! Author: Frederic Hourdin, January 30th, 1992
12
13	! Object: computation of the soil temperature evolution, the heat
14	! capacity per unit surface and the surface conduction flux
15
16	! Method: implicit time integration
17
18	! Consecutive ground temperatures are related by:
19	! 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	! 1) new temperatures are computed using (1)
23	! 2) C and D coefficients are computed from the new temperature
24	! profile for the t + dt time-step
25	! 3) the coefficients A and B are computed where the diffusive
26	! 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
34	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
39	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
44	real, intent(inout):: tsoil(:, :) ! (knon, nsoilmx)
45	! temperature inside the ground (K)
46
47	REAL, intent(out):: soilcap(:) ! (knon)
48	! specific heat per unit surface (W m-2 s K-1)
49
50	REAL, intent(out):: soilflux(:) ! (knon)
51	! surface diffusive flux from ground (W m-2)
52
53	! Local:
54
55	INTEGER knon, ig, jk
56	REAL zdz2(nsoilmx)
57	real z1(size(tsurf), nbsrf) ! (knon, nbsrf)
58	REAL min_period, dalph_soil
59	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
66	! Depths:
67	REAL rk, fz1, rk1, rk2
68
69	!-----------------------------------------------------------------------
70
71	knon = size(tsurf)
72
73	! 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
78	IF (nisurf==is_sic) THEN
79	DO ig = 1, knon
80	ztherm_i(ig) = iice
81	IF (snow(ig) > 0.0) ztherm_i(ig) = isno
82	END DO
83	ELSE IF (nisurf==is_lic) THEN
84	DO ig = 1, knon
85	ztherm_i(ig) = iice
86	IF (snow(ig) > 0.0) ztherm_i(ig) = isno
87	END DO
88	ELSE IF (nisurf==is_ter) THEN
89	DO ig = 1, knon
90	ztherm_i(ig) = isol
91	IF (snow(ig) > 0.0) ztherm_i(ig) = isno
92	END DO
93	ELSE IF (nisurf==is_oce) THEN
94	DO ig = 1, knon
95	ztherm_i(ig) = iice
96	END DO
97	ELSE
98	PRINT *, 'valeur d indice non prevue', nisurf
99	STOP 1
100	END IF
101
102	IF (firstsurf(nisurf)) THEN
103	! ground levels
104	! grnd=z / l where l is the skin depth of the diurnal cycle:
105
106	min_period = 1800. ! en secondes
107	dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
108
109	OPEN(99, FILE='soil.def', STATUS='old', FORM='formatted', ERR=9999)
110	READ(99, *) min_period
111	READ(99, *) dalph_soil
112	PRINT *, 'Discretization for the soil model'
113	PRINT *, 'First level e-folding depth', min_period, ' dalph', &
114	dalph_soil
115	CLOSE(99)
116	9999 CONTINUE
117
118	! la premiere couche represente un dixieme de cycle diurne
119	fz1 = sqrt(min_period / 3.14)
120
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	dz1(jk) = 1. / (fz(rk1) - fz(rk2))
130	END DO
131	lambda = fz(.5) * dz1(1)
132	PRINT *, 'full layers, intermediate layers (seconds)'
133	DO jk = 1, nsoilmx
134	rk = jk
135	rk1 = jk + .5
136	rk2 = jk - .5
137	PRINT , 'fz=', fz(rk1) fz(rk2) * 3.14, fz(rk) * fz(rk) * 3.14
138	END DO
139	! PB
140	firstsurf(nisurf) = .FALSE.
141	ELSE
142	! Computation of the soil temperatures using the Zc and Zd
143	! coefficient computed at the previous time-step:
144
145	! surface temperature
146	DO ig = 1, knon
147	tsoil(ig, 1) = (lambda * zc(ig, 1, nisurf) + tsurf(ig)) &
148	/ (lambda * (1. - zd(ig, 1, nisurf)) + 1.)
149	END DO
150
151	! other temperatures
152	DO jk = 1, nsoilmx - 1
153	DO ig = 1, knon
154	tsoil(ig, jk + 1) = zc(ig, jk, nisurf) &
155	+ zd(ig, jk, nisurf) * tsoil(ig, jk)
156	END DO
157	END DO
158	END IF
159
160	! Computation of the Zc and Zd coefficient for the next step:
161
162	IF (nisurf==is_sic) THEN
163	DO ig = 1, knon
164	tsoil(ig, nsoilmx) = rtt - 1.8
165	END DO
166	END IF
167
168	DO jk = 1, nsoilmx
169	zdz2(jk) = dz2(jk) / dtime
170	END DO
171
172	DO ig = 1, knon
173	z1(ig, nisurf) = zdz2(nsoilmx) + dz1(nsoilmx - 1)
174	zc(ig, nsoilmx - 1, nisurf) = zdz2(nsoilmx) * tsoil(ig, nsoilmx) / &
175	z1(ig, nisurf)
176	zd(ig, nsoilmx - 1, nisurf) = dz1(nsoilmx - 1) / z1(ig, nisurf)
177	END DO
178
179	DO jk = nsoilmx - 1, 2, - 1
180	DO ig = 1, knon
181	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	END DO
187	END DO
188
189	! computation of the surface diffusive flux from ground and
190	! calorific capacity of the ground:
191
192	DO ig = 1, knon
193	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	END DO
202
203	contains
204
205	pure real function fz(rk)
206
207	real, intent(in):: rk
208
209	!-----------------------------------------
210
211	fz = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
212
213	end function fz
214
215	END SUBROUTINE soil
216
217	end module soil_m