phylmd/Interface_surf/soil.f

module soil_m

  IMPLICIT NONE

contains

  SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil, pcapcal, &
       pfluxgrd)

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

    USE dimens_m
    USE indicesol
    USE dimphy
    USE dimsoil
    USE suphec_m

    ! =======================================================================

    ! Auteur:  Frederic Hourdin     30/01/92
    ! -------

    ! objet:  computation of : the soil temperature evolution
    ! ------                   the surfacic heat capacity "Capcal"
    ! the surface conduction flux pcapcal


    ! 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 + Capcal (Ts(t+dt)-Ts(t))/dt
    ! with F0 = A + B (Ts(t))
    ! Capcal = B*dt

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

    ! Arguments:
    ! ----------
    ! ptimestep            physical timestep (s)
    ! indice               sub-surface index
    ! snow(klon,nbsrf)     snow
    ! ptsrf(klon)          surface temperature at time-step t (K)
    ! ptsoil(klon,nsoilmx) temperature inside the ground (K)
    ! pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
    ! pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)

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


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

    REAL ptimestep
    INTEGER indice, knon
    REAL ptsrf(klon), ptsoil(klon, nsoilmx), snow(klon)
    REAL pcapcal(klon), pfluxgrd(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 firstcall, firstsurf(nbsrf)
    SAVE firstcall, firstsurf
    REAL isol, isno, iice
    SAVE isol, isno, iice

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

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

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

    REAL fz, rk, fz1, rk1, rk2

    fz(rk) = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
    pfluxgrd(:) = 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 (indice==is_sic) THEN
       DO ig = 1, knon
          ztherm_i(ig) = iice
          IF (snow(ig)>0.0) ztherm_i(ig) = isno
       END DO
    ELSE IF (indice==is_lic) THEN
       DO ig = 1, knon
          ztherm_i(ig) = iice
          IF (snow(ig)>0.0) ztherm_i(ig) = isno
       END DO
    ELSE IF (indice==is_ter) THEN
       DO ig = 1, knon
          ztherm_i(ig) = isol
          IF (snow(ig)>0.0) ztherm_i(ig) = isno
       END DO
    ELSE IF (indice==is_oce) THEN
       DO ig = 1, knon
          ztherm_i(ig) = iice
       END DO
    ELSE
       PRINT *, 'valeur d indice non prevue', indice
       STOP 1
    END IF


    ! $$$      IF (firstcall) THEN
    IF (firstsurf(indice)) 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(indice) = .FALSE.
       ! $$$         firstcall =.false.

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

    ELSE !--not firstcall
       ! -----------------------------------------------------------------------
       ! Computation of the soil temperatures using the Cgrd and Dgrd
       ! coefficient computed at the previous time-step:
       ! -----------------------------------------------

       ! surface temperature
       DO ig = 1, knon
          ptsoil(ig, 1) = (lambda*zc(ig,1,indice)+ptsrf(ig))/(lambda*(1.-zd(ig,1, &
               indice))+1.)
       END DO

       ! other temperatures
       DO jk = 1, nsoilmx - 1
          DO ig = 1, knon
             ptsoil(ig, jk+1) = zc(ig, jk, indice) + zd(ig, jk, indice)*ptsoil(ig, &
                  jk)
          END DO
       END DO

    END IF !--not firstcall
    ! -----------------------------------------------------------------------
    ! Computation of the Cgrd and Dgrd coefficient for the next step:
    ! ---------------------------------------------------------------

    ! $$$  PB ajout pour cas glace de mer
    IF (indice==is_sic) THEN
       DO ig = 1, knon
          ptsoil(ig, nsoilmx) = rtt - 1.8
       END DO
    END IF

    DO jk = 1, nsoilmx
       zdz2(jk) = dz2(jk)/ptimestep
    END DO

    DO ig = 1, knon
       z1(ig, indice) = zdz2(nsoilmx) + dz1(nsoilmx-1)
       zc(ig, nsoilmx-1, indice) = zdz2(nsoilmx)*ptsoil(ig, nsoilmx)/ &
            z1(ig, indice)
       zd(ig, nsoilmx-1, indice) = dz1(nsoilmx-1)/z1(ig, indice)
    END DO

    DO jk = nsoilmx - 1, 2, -1
       DO ig = 1, knon
          z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice)))
          zc(ig, jk-1, indice) = (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice) &
               )*z1(ig, indice)
          zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice)
       END DO
    END DO

    ! -----------------------------------------------------------------------
    ! computation of the surface diffusive flux from ground and
    ! calorific capacity of the ground:
    ! ---------------------------------

    DO ig = 1, knon
       pfluxgrd(ig) = ztherm_i(ig)*dz1(1)*(zc(ig,1,indice)+(zd(ig,1, &
            indice)-1.)*ptsoil(ig,1))
       pcapcal(ig) = ztherm_i(ig)*(dz2(1)+ptimestep*(1.-zd(ig,1,indice))*dz1(1))
       z1(ig, indice) = lambda*(1.-zd(ig,1,indice)) + 1.
       pcapcal(ig) = pcapcal(ig)/z1(ig, indice)
       pfluxgrd(ig) = pfluxgrd(ig) + pcapcal(ig)*(ptsoil(ig,1)*z1(ig,indice)- &
            lambda*zc(ig,1,indice)-ptsrf(ig))/ptimestep
    END DO

  END SUBROUTINE soil

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