phylmd/Interface_surf/soil.f90

module soil_m

  IMPLICIT NONE

  private fz

contains

  SUBROUTINE soil(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) (equation 1)
    ! The coefficients C and D are computed at the t - dt time-step.
    ! Structure of the procedure:
    ! 1) new temperatures are computed using equation 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

    ! Libraries:
    use jumble, only: new_unit

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

    INTEGER, intent(in):: nisurf ! surface type 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), layer 1 nearest to the surface

    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, unit
    REAL zdz2(nsoilmx)
    real z1(size(tsurf), nbsrf) ! (knon, nbsrf)
    REAL min_period ! in s
    real dalph_soil ! rapport entre les \'epaisseurs de 2 couches successives
    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, parameter:: isol = 2000., isno = 2000., iice = 2000.
    REAL rk, fz1, rk1, rk2 ! depths

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

    knon = size(tsurf)

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

    select case (nisurf)
    case (is_sic)
       DO ig = 1, knon
          IF (snow(ig) > 0.) then
             ztherm_i(ig) = isno
          else
             ztherm_i(ig) = iice
          end IF
       END DO
    case (is_lic)
       DO ig = 1, knon
          IF (snow(ig) > 0.) then
             ztherm_i(ig) = isno
          else
             ztherm_i(ig) = iice
          end IF
       END DO
    case (is_ter)
       DO ig = 1, knon
          IF (snow(ig) > 0.) then
             ztherm_i(ig) = isno
          else
             ztherm_i(ig) = isol
          end IF
       END DO
    case (is_oce)
       DO ig = 1, knon
          ztherm_i(ig) = iice
       END DO
    case default
       PRINT *, 'soil: unexpected subscript value:', nisurf
       STOP 1
    END select

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

       min_period = 1800.
       dalph_soil = 2.
       call new_unit(unit)
       OPEN(unit, FILE = 'soil.def', STATUS = 'old', action = "read", &
            position = 'rewind', ERR = 9999)
       READ(unit, fmt = *) min_period
       READ(unit, fmt = *) dalph_soil
       PRINT *, 'Discretization for the soil model'
       PRINT *, 'First level e-folding depth', min_period, ' dalph', &
            dalph_soil
       CLOSE(unit)
9999   CONTINUE

       ! La premi\`ere couche repr\'esente un dixi\`eme de cycle diurne :
       fz1 = sqrt(min_period / 3.14)

       DO jk = 1, nsoilmx
          rk1 = jk
          rk2 = jk - 1
          dz2(jk) = fz(rk1, dalph_soil, fz1) - fz(rk2, dalph_soil, fz1)
       END DO

       DO jk = 1, nsoilmx - 1
          rk1 = jk + .5
          rk2 = jk - .5
          dz1(jk) = 1. / (fz(rk1, dalph_soil, fz1) - fz(rk2, dalph_soil, fz1))
       END DO

       lambda = fz(.5, dalph_soil, fz1) * dz1(1)
       PRINT *, 'full layers, intermediate layers (seconds)'

       DO jk = 1, nsoilmx
          rk = jk
          rk1 = jk + .5
          rk2 = jk - .5
          PRINT *, 'fz=', fz(rk1, dalph_soil, fz1) * fz(rk2, dalph_soil, fz1) &
               * 3.14, fz(rk, dalph_soil, fz1) * fz(rk, dalph_soil, fz1) * 3.14
       END DO

       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) / dtphys
    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) &
            + dtphys * (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)) / dtphys
    END DO

  END SUBROUTINE soil

  !****************************************************************

  pure real function fz(rk, dalph_soil, fz1)

    real, intent(in):: rk

    real, intent(in):: dalph_soil
    ! rapport entre les \'epaisseurs de 2 couches successives

    real, intent(in):: fz1 ! depth

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

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

  end function fz

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