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