phylmd/Interface_surf/soil.f

module soil_m

  IMPLICIT NONE

contains

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

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

    ! Author: Frederic Hourdin 30/01/92

    ! Object: computation of the soil temperature evolution, the
    ! surfacic heat capacity "Soilcap" 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.
    ! Routine structure:
    ! 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 dimens_m, only:
    USE indicesol
    USE dimphy
    USE dimsoil
    USE suphec_m

    ! Interface:
    ! ----------

    ! Arguments:
    ! ----------
    ! dtime physical timestep (s)
    ! indice sub-surface index
    ! snow(klon, nbsrf) snow
    ! tsurf(knon) surface temperature at time-step t (K)
    ! tsoil(klon, nsoilmx) temperature inside the ground (K)
    ! soilcap(klon) surfacic specific heat (W*m-2*s*K-1)
    ! soilflux(klon) surface diffusive flux from ground (Wm-2)

    ! declarations:
    ! -------------

    ! -----------------------------------------------------------------------
    ! arguments
    ! ---------

    REAL dtime
    INTEGER nisurf, knon
    REAL tsurf(knon), tsoil(klon, nsoilmx), snow(klon)
    REAL soilcap(klon), soilflux(klon)

    ! -----------------------------------------------------------------------
    ! local arrays
    ! ------------

    INTEGER ig, jk
    ! $$$ REAL zdz2(nsoilmx), z1(klon)
    REAL zdz2(nsoilmx), z1(klon, nbsrf)
    REAL min_period, dalph_soil
    REAL ztherm_i(klon)

    ! local saved variables:
    ! ----------------------
    REAL dz1(nsoilmx), dz2(nsoilmx)
    ! $$$ REAL zc(klon, nsoilmx), zd(klon, nsoilmx)
    REAL zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf)
    REAL lambda
    SAVE dz1, dz2, zc, zd, lambda
    LOGICAL firstsurf(nbsrf)
    SAVE firstsurf
    REAL isol, isno, iice
    SAVE isol, isno, iice

    DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./

    DATA isol, isno, iice/2000., 2000., 2000./

    ! -----------------------------------------------------------------------
    ! Depthts:
    ! --------

    REAL rk, fz1, rk1, rk2

    soilflux(:) = 0.
    ! calcul de l'inertie thermique a partir de la variable rnat.
    ! on initialise a iice meme au-dessus d'un point de mer au cas
    ! ou 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.

       ! Initialisations:
       ! ----------------

    ELSE
       ! -----------------------------------------------------------------------
       ! Computation of the soil temperatures using the Cgrd and Dgrd
       ! 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 Cgrd and Dgrd coefficient for the next step:
    ! ---------------------------------------------------------------

    ! $$$ PB ajout pour cas glace de mer
    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

    real function fz(rk)
      real 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, knon, 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 30/01/92
12
13	! Object: computation of the soil temperature evolution, the
14	! surfacic heat capacity "Soilcap" 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	! Routine structure:
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 Fdiff = F0 + Soilcap (Ts(t+dt)-Ts(t))/dt
29	! with F0 = A + B (Ts(t))
30	! Soilcap = B*dt
31
32	USE dimens_m, only:
33	USE indicesol
34	USE dimphy
35	USE dimsoil
36	USE suphec_m
37
38	! Interface:
39	! ----------
40
41	! Arguments:
42	! ----------
43	! dtime physical timestep (s)
44	! indice sub-surface index
45	! snow(klon, nbsrf) snow
46	! tsurf(knon) surface temperature at time-step t (K)
47	! tsoil(klon, nsoilmx) temperature inside the ground (K)
48	! soilcap(klon) surfacic specific heat (Wm-2s*K-1)
49	! soilflux(klon) surface diffusive flux from ground (Wm-2)
50
51	! declarations:
52	! -------------
53
54	! -----------------------------------------------------------------------
55	! arguments
56	! ---------
57
58	REAL dtime
59	INTEGER nisurf, knon
60	REAL tsurf(knon), tsoil(klon, nsoilmx), snow(klon)
61	REAL soilcap(klon), soilflux(klon)
62
63	! -----------------------------------------------------------------------
64	! local arrays
65	! ------------
66
67	INTEGER ig, jk
68	! $$$ REAL zdz2(nsoilmx), z1(klon)
69	REAL zdz2(nsoilmx), z1(klon, nbsrf)
70	REAL min_period, dalph_soil
71	REAL ztherm_i(klon)
72
73	! local saved variables:
74	! ----------------------
75	REAL dz1(nsoilmx), dz2(nsoilmx)
76	! $$$ REAL zc(klon, nsoilmx), zd(klon, nsoilmx)
77	REAL zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf)
78	REAL lambda
79	SAVE dz1, dz2, zc, zd, lambda
80	LOGICAL firstsurf(nbsrf)
81	SAVE firstsurf
82	REAL isol, isno, iice
83	SAVE isol, isno, iice
84
85	DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./
86
87	DATA isol, isno, iice/2000., 2000., 2000./
88
89	! -----------------------------------------------------------------------
90	! Depthts:
91	! --------
92
93	REAL rk, fz1, rk1, rk2
94
95	soilflux(:) = 0.
96	! calcul de l'inertie thermique a partir de la variable rnat.
97	! on initialise a iice meme au-dessus d'un point de mer au cas
98	! ou le point de mer devienne point de glace au pas suivant
99	! on corrige si on a un point de terre avec ou sans glace
100
101	IF (nisurf==is_sic) THEN
102	DO ig = 1, knon
103	ztherm_i(ig) = iice
104	IF (snow(ig)>0.0) ztherm_i(ig) = isno
105	END DO
106	ELSE IF (nisurf==is_lic) THEN
107	DO ig = 1, knon
108	ztherm_i(ig) = iice
109	IF (snow(ig)>0.0) ztherm_i(ig) = isno
110	END DO
111	ELSE IF (nisurf==is_ter) THEN
112	DO ig = 1, knon
113	ztherm_i(ig) = isol
114	IF (snow(ig)>0.0) ztherm_i(ig) = isno
115	END DO
116	ELSE IF (nisurf==is_oce) THEN
117	DO ig = 1, knon
118	ztherm_i(ig) = iice
119	END DO
120	ELSE
121	PRINT *, 'valeur d indice non prevue', nisurf
122	STOP 1
123	END IF
124
125	IF (firstsurf(nisurf)) THEN
126
127	! -----------------------------------------------------------------------
128	! ground levels
129	! grnd=z/l where l is the skin depth of the diurnal cycle:
130	! --------------------------------------------------------
131
132	min_period = 1800. ! en secondes
133	dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
134
135	OPEN (99, FILE='soil.def', STATUS='old', FORM='formatted', ERR=9999)
136	READ (99, *) min_period
137	READ (99, *) dalph_soil
138	PRINT *, 'Discretization for the soil model'
139	PRINT *, 'First level e-folding depth', min_period, ' dalph', &
140	dalph_soil
141	CLOSE (99)
142	9999 CONTINUE
143
144	! la premiere couche represente un dixieme de cycle diurne
145	fz1 = sqrt(min_period/3.14)
146
147	DO jk = 1, nsoilmx
148	rk1 = jk
149	rk2 = jk - 1
150	dz2(jk) = fz(rk1) - fz(rk2)
151	END DO
152	DO jk = 1, nsoilmx - 1
153	rk1 = jk + .5
154	rk2 = jk - .5
155	dz1(jk) = 1./(fz(rk1)-fz(rk2))
156	END DO
157	lambda = fz(.5)*dz1(1)
158	PRINT *, 'full layers, intermediate layers (seconds)'
159	DO jk = 1, nsoilmx
160	rk = jk
161	rk1 = jk + .5
162	rk2 = jk - .5
163	PRINT , 'fz=', fz(rk1)fz(rk2)3.14, fz(rk)fz(rk)*3.14
164	END DO
165	! PB
166	firstsurf(nisurf) = .FALSE.
167
168	! Initialisations:
169	! ----------------
170
171	ELSE
172	! -----------------------------------------------------------------------
173	! Computation of the soil temperatures using the Cgrd and Dgrd
174	! coefficient computed at the previous time-step:
175	! -----------------------------------------------
176
177	! surface temperature
178	DO ig = 1, knon
179	tsoil(ig, 1) = (lambdazc(ig, 1, nisurf)+tsurf(ig))/(lambda(1.-zd(ig, 1, &
180	nisurf))+1.)
181	END DO
182
183	! other temperatures
184	DO jk = 1, nsoilmx - 1
185	DO ig = 1, knon
186	tsoil(ig, jk+1) = zc(ig, jk, nisurf) + zd(ig, jk, nisurf)*tsoil(ig, &
187	jk)
188	END DO
189	END DO
190
191	END IF
192	! -----------------------------------------------------------------------
193	! Computation of the Cgrd and Dgrd coefficient for the next step:
194	! ---------------------------------------------------------------
195
196	! $$$ PB ajout pour cas glace de mer
197	IF (nisurf==is_sic) THEN
198	DO ig = 1, knon
199	tsoil(ig, nsoilmx) = rtt - 1.8
200	END DO
201	END IF
202
203	DO jk = 1, nsoilmx
204	zdz2(jk) = dz2(jk)/dtime
205	END DO
206
207	DO ig = 1, knon
208	z1(ig, nisurf) = zdz2(nsoilmx) + dz1(nsoilmx-1)
209	zc(ig, nsoilmx-1, nisurf) = zdz2(nsoilmx)*tsoil(ig, nsoilmx)/ &
210	z1(ig, nisurf)
211	zd(ig, nsoilmx-1, nisurf) = dz1(nsoilmx-1)/z1(ig, nisurf)
212	END DO
213
214	DO jk = nsoilmx - 1, 2, -1
215	DO ig = 1, knon
216	z1(ig, nisurf) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig, jk, nisurf)))
217	zc(ig, jk-1, nisurf) = (tsoil(ig, jk)zdz2(jk)+dz1(jk)zc(ig, jk, nisurf) &
218	)*z1(ig, nisurf)
219	zd(ig, jk-1, nisurf) = dz1(jk-1)*z1(ig, nisurf)
220	END DO
221	END DO
222
223	! -----------------------------------------------------------------------
224	! computation of the surface diffusive flux from ground and
225	! calorific capacity of the ground:
226	! ---------------------------------
227
228	DO ig = 1, knon
229	soilflux(ig) = ztherm_i(ig)dz1(1)(zc(ig, 1, nisurf)+(zd(ig, 1, &
230	nisurf)-1.)*tsoil(ig, 1))
231	soilcap(ig) = ztherm_i(ig)(dz2(1)+dtime(1.-zd(ig, 1, nisurf))*dz1(1))
232	z1(ig, nisurf) = lambda*(1.-zd(ig, 1, nisurf)) + 1.
233	soilcap(ig) = soilcap(ig)/z1(ig, nisurf)
234	soilflux(ig) = soilflux(ig) + soilcap(ig)(tsoil(ig, 1)z1(ig, nisurf)- &
235	lambda*zc(ig, 1, nisurf)-tsurf(ig))/dtime
236	END DO
237
238	contains
239
240	real function fz(rk)
241	real rk
242	fz = fz1(dalph_soil*rk-1.)/(dalph_soil-1.)
243	end function fz
244
245	END SUBROUTINE soil
246
247	end module soil_m