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	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 indicesol
33			USE dimphy
34			USE dimsoil
35			USE suphec_m
36
37	guez	101	! Interface:
38			! ----------
39	guez	3
40	guez	101	! Arguments:
41			! ----------
42	guez	175	! 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	guez	3
50	guez	101	! declarations:
51			! -------------
52	guez	3
53	guez	101	! -----------------------------------------------------------------------
54			! arguments
55			! ---------
56	guez	3
57	guez	175	REAL dtime
58			INTEGER nisurf, knon
59			REAL tsurf(knon), tsoil(klon, nsoilmx), snow(klon)
60			REAL soilcap(klon), soilflux(klon)
61	guez	3
62	guez	101	! -----------------------------------------------------------------------
63			! local arrays
64			! ------------
65	guez	3
66	guez	101	INTEGER ig, jk
67	guez	175	! $$$ REAL zdz2(nsoilmx), z1(klon)
68	guez	101	REAL zdz2(nsoilmx), z1(klon, nbsrf)
69			REAL min_period, dalph_soil
70			REAL ztherm_i(klon)
71	guez	3
72	guez	101	! local saved variables:
73			! ----------------------
74			REAL dz1(nsoilmx), dz2(nsoilmx)
75	guez	175	! $$$ REAL zc(klon, nsoilmx), zd(klon, nsoilmx)
76	guez	101	REAL zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf)
77			REAL lambda
78			SAVE dz1, dz2, zc, zd, lambda
79	guez	157	LOGICAL firstsurf(nbsrf)
80			SAVE firstsurf
81	guez	101	REAL isol, isno, iice
82			SAVE isol, isno, iice
83	guez	3
84	guez	101	DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./
85	guez	3
86	guez	101	DATA isol, isno, iice/2000., 2000., 2000./
87	guez	3
88	guez	101	! -----------------------------------------------------------------------
89			! Depthts:
90			! --------
91	guez	3
92	guez	108	REAL rk, fz1, rk1, rk2
93	guez	3
94	guez	175	soilflux(:) = 0.
95	guez	101	! 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	guez	3
100	guez	175	IF (nisurf==is_sic) THEN
101	guez	101	DO ig = 1, knon
102			ztherm_i(ig) = iice
103			IF (snow(ig)>0.0) ztherm_i(ig) = isno
104			END DO
105	guez	175	ELSE IF (nisurf==is_lic) THEN
106	guez	101	DO ig = 1, knon
107			ztherm_i(ig) = iice
108			IF (snow(ig)>0.0) ztherm_i(ig) = isno
109			END DO
110	guez	175	ELSE IF (nisurf==is_ter) THEN
111	guez	101	DO ig = 1, knon
112			ztherm_i(ig) = isol
113			IF (snow(ig)>0.0) ztherm_i(ig) = isno
114			END DO
115	guez	175	ELSE IF (nisurf==is_oce) THEN
116	guez	101	DO ig = 1, knon
117			ztherm_i(ig) = iice
118			END DO
119			ELSE
120	guez	175	PRINT *, 'valeur d indice non prevue', nisurf
121	guez	101	STOP 1
122			END IF
123	guez	3
124	guez	175	IF (firstsurf(nisurf)) THEN
125	guez	3
126	guez	101	! -----------------------------------------------------------------------
127			! ground levels
128			! grnd=z/l where l is the skin depth of the diurnal cycle:
129			! --------------------------------------------------------
130	guez	3
131	guez	101	min_period = 1800. ! en secondes
132			dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
133	guez	81
134	guez	101	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	guez	175	PRINT *, 'First level e-folding depth', min_period, ' dalph', &
139	guez	101	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	guez	175	firstsurf(nisurf) = .FALSE.
166	guez	101
167			! Initialisations:
168			! ----------------
169
170	guez	157	ELSE
171	guez	101	! -----------------------------------------------------------------------
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	guez	175	tsoil(ig, 1) = (lambdazc(ig, 1, nisurf)+tsurf(ig))/(lambda(1.-zd(ig, 1, &
179			nisurf))+1.)
180	guez	101	END DO
181
182			! other temperatures
183			DO jk = 1, nsoilmx - 1
184			DO ig = 1, knon
185	guez	175	tsoil(ig, jk+1) = zc(ig, jk, nisurf) + zd(ig, jk, nisurf)*tsoil(ig, &
186	guez	101	jk)
187			END DO
188			END DO
189
190	guez	157	END IF
191	guez	81	! -----------------------------------------------------------------------
192	guez	101	! Computation of the Cgrd and Dgrd coefficient for the next step:
193			! ---------------------------------------------------------------
194	guez	81
195	guez	175	! $$$ PB ajout pour cas glace de mer
196			IF (nisurf==is_sic) THEN
197	guez	101	DO ig = 1, knon
198	guez	175	tsoil(ig, nsoilmx) = rtt - 1.8
199	guez	101	END DO
200			END IF
201	guez	81
202	guez	101	DO jk = 1, nsoilmx
203	guez	175	zdz2(jk) = dz2(jk)/dtime
204	guez	81	END DO
205
206			DO ig = 1, knon
207	guez	175	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	guez	81	END DO
212
213	guez	101	DO jk = nsoilmx - 1, 2, -1
214			DO ig = 1, knon
215	guez	175	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	guez	101	END DO
220			END DO
221	guez	81
222	guez	101	! -----------------------------------------------------------------------
223			! computation of the surface diffusive flux from ground and
224			! calorific capacity of the ground:
225			! ---------------------------------
226	guez	81
227			DO ig = 1, knon
228	guez	175	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	guez	81	END DO
236
237	guez	108	contains
238
239			real function fz(rk)
240	guez	202
241			real, intent(in):: rk
242
243			!-----------------------------------------
244
245	guez	108	fz = fz1(dalph_soil*rk-1.)/(dalph_soil-1.)
246	guez	202
247	guez	108	end function fz
248
249	guez	101	END SUBROUTINE soil
250	guez	81
251	guez	101	end module soil_m