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