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	guez	101	module soil_m
2	guez	3
3	guez	81	IMPLICIT NONE
4	guez	3
5	guez	101	contains
6	guez	3
7	guez	175	SUBROUTINE soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, soilflux)
8	guez	3
9	guez	101	! From LMDZ4/libf/phylmd/soil.F, version 1.1.1.1 2004/05/19
10	guez	3
11	guez	175	! Author: Frederic Hourdin 30/01/92
12	guez	3
13	guez	175	! Object: computation of the soil temperature evolution, the
14			! surfacic heat capacity "Soilcap" and the surface conduction flux
15	guez	3
16	guez	175	! Method: implicit time integration
17	guez	3
18	guez	101	! Consecutive ground temperatures are related by:
19	guez	175	! T(k+1) = C(k) + D(k)*T(k) (1)
20	guez	101	! the coefficients C and D are computed at the t-dt time-step.
21			! Routine structure:
22	guez	175	! 1) new temperatures are computed using (1)
23			! 2) C and D coefficients are computed from the new temperature
24	guez	101	! profile for the t+dt time-step
25	guez	175	! 3) the coefficients A and B are computed where the diffusive
26	guez	101	! fluxes at the t+dt time-step is given by
27			! Fdiff = A + B Ts(t+dt)
28	guez	175	! or Fdiff = F0 + Soilcap (Ts(t+dt)-Ts(t))/dt
29	guez	101	! with F0 = A + B (Ts(t))
30	guez	175	! Soilcap = B*dt
31	guez	3
32	guez	175	USE dimens_m, only:
33			USE indicesol
34			USE dimphy
35			USE dimsoil
36			USE suphec_m
37
38	guez	101	! Interface:
39			! ----------
40	guez	3
41	guez	101	! Arguments:
42			! ----------
43	guez	175	! 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	guez	3
51	guez	101	! declarations:
52			! -------------
53	guez	3
54	guez	101	! -----------------------------------------------------------------------
55			! arguments
56			! ---------
57	guez	3
58	guez	175	REAL dtime
59			INTEGER nisurf, knon
60			REAL tsurf(knon), tsoil(klon, nsoilmx), snow(klon)
61			REAL soilcap(klon), soilflux(klon)
62	guez	3
63	guez	101	! -----------------------------------------------------------------------
64			! local arrays
65			! ------------
66	guez	3
67	guez	101	INTEGER ig, jk
68	guez	175	! $$$ REAL zdz2(nsoilmx), z1(klon)
69	guez	101	REAL zdz2(nsoilmx), z1(klon, nbsrf)
70			REAL min_period, dalph_soil
71			REAL ztherm_i(klon)
72	guez	3
73	guez	101	! local saved variables:
74			! ----------------------
75			REAL dz1(nsoilmx), dz2(nsoilmx)
76	guez	175	! $$$ REAL zc(klon, nsoilmx), zd(klon, nsoilmx)
77	guez	101	REAL zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf)
78			REAL lambda
79			SAVE dz1, dz2, zc, zd, lambda
80	guez	157	LOGICAL firstsurf(nbsrf)
81			SAVE firstsurf
82	guez	101	REAL isol, isno, iice
83			SAVE isol, isno, iice
84	guez	3
85	guez	101	DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./
86	guez	3
87	guez	101	DATA isol, isno, iice/2000., 2000., 2000./
88	guez	3
89	guez	101	! -----------------------------------------------------------------------
90			! Depthts:
91			! --------
92	guez	3
93	guez	108	REAL rk, fz1, rk1, rk2
94	guez	3
95	guez	175	soilflux(:) = 0.
96	guez	101	! 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	guez	3
101	guez	175	IF (nisurf==is_sic) THEN
102	guez	101	DO ig = 1, knon
103			ztherm_i(ig) = iice
104			IF (snow(ig)>0.0) ztherm_i(ig) = isno
105			END DO
106	guez	175	ELSE IF (nisurf==is_lic) THEN
107	guez	101	DO ig = 1, knon
108			ztherm_i(ig) = iice
109			IF (snow(ig)>0.0) ztherm_i(ig) = isno
110			END DO
111	guez	175	ELSE IF (nisurf==is_ter) THEN
112	guez	101	DO ig = 1, knon
113			ztherm_i(ig) = isol
114			IF (snow(ig)>0.0) ztherm_i(ig) = isno
115			END DO
116	guez	175	ELSE IF (nisurf==is_oce) THEN
117	guez	101	DO ig = 1, knon
118			ztherm_i(ig) = iice
119			END DO
120			ELSE
121	guez	175	PRINT *, 'valeur d indice non prevue', nisurf
122	guez	101	STOP 1
123			END IF
124	guez	3
125	guez	175	IF (firstsurf(nisurf)) THEN
126	guez	3
127	guez	101	! -----------------------------------------------------------------------
128			! ground levels
129			! grnd=z/l where l is the skin depth of the diurnal cycle:
130			! --------------------------------------------------------
131	guez	3
132	guez	101	min_period = 1800. ! en secondes
133			dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
134	guez	81
135	guez	101	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	guez	175	PRINT *, 'First level e-folding depth', min_period, ' dalph', &
140	guez	101	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	guez	175	firstsurf(nisurf) = .FALSE.
167	guez	101
168			! Initialisations:
169			! ----------------
170
171	guez	157	ELSE
172	guez	101	! -----------------------------------------------------------------------
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	guez	175	tsoil(ig, 1) = (lambdazc(ig, 1, nisurf)+tsurf(ig))/(lambda(1.-zd(ig, 1, &
180			nisurf))+1.)
181	guez	101	END DO
182
183			! other temperatures
184			DO jk = 1, nsoilmx - 1
185			DO ig = 1, knon
186	guez	175	tsoil(ig, jk+1) = zc(ig, jk, nisurf) + zd(ig, jk, nisurf)*tsoil(ig, &
187	guez	101	jk)
188			END DO
189			END DO
190
191	guez	157	END IF
192	guez	81	! -----------------------------------------------------------------------
193	guez	101	! Computation of the Cgrd and Dgrd coefficient for the next step:
194			! ---------------------------------------------------------------
195	guez	81
196	guez	175	! $$$ PB ajout pour cas glace de mer
197			IF (nisurf==is_sic) THEN
198	guez	101	DO ig = 1, knon
199	guez	175	tsoil(ig, nsoilmx) = rtt - 1.8
200	guez	101	END DO
201			END IF
202	guez	81
203	guez	101	DO jk = 1, nsoilmx
204	guez	175	zdz2(jk) = dz2(jk)/dtime
205	guez	81	END DO
206
207			DO ig = 1, knon
208	guez	175	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	guez	81	END DO
213
214	guez	101	DO jk = nsoilmx - 1, 2, -1
215			DO ig = 1, knon
216	guez	175	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	guez	101	END DO
221			END DO
222	guez	81
223	guez	101	! -----------------------------------------------------------------------
224			! computation of the surface diffusive flux from ground and
225			! calorific capacity of the ground:
226			! ---------------------------------
227	guez	81
228			DO ig = 1, knon
229	guez	175	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	guez	81	END DO
237
238	guez	108	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	guez	101	END SUBROUTINE soil
246	guez	81
247	guez	101	end module soil_m