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(knon)          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(knon), 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 rk, fz1, rk1, rk2

    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

  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(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(knon) 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(knon), 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 rk, fz1, rk1, rk2
103
104	pfluxgrd(:) = 0.
105	! calcul de l'inertie thermique a partir de la variable rnat.
106	! on initialise a iice meme au-dessus d'un point de mer au cas
107	! ou le point de mer devienne point de glace au pas suivant
108	! on corrige si on a un point de terre avec ou sans glace
109
110	IF (indice==is_sic) THEN
111	DO ig = 1, knon
112	ztherm_i(ig) = iice
113	IF (snow(ig)>0.0) ztherm_i(ig) = isno
114	END DO
115	ELSE IF (indice==is_lic) THEN
116	DO ig = 1, knon
117	ztherm_i(ig) = iice
118	IF (snow(ig)>0.0) ztherm_i(ig) = isno
119	END DO
120	ELSE IF (indice==is_ter) THEN
121	DO ig = 1, knon
122	ztherm_i(ig) = isol
123	IF (snow(ig)>0.0) ztherm_i(ig) = isno
124	END DO
125	ELSE IF (indice==is_oce) THEN
126	DO ig = 1, knon
127	ztherm_i(ig) = iice
128	END DO
129	ELSE
130	PRINT *, 'valeur d indice non prevue', indice
131	STOP 1
132	END IF
133
134
135	! $$$ IF (firstcall) THEN
136	IF (firstsurf(indice)) THEN
137
138	! -----------------------------------------------------------------------
139	! ground levels
140	! grnd=z/l where l is the skin depth of the diurnal cycle:
141	! --------------------------------------------------------
142
143	min_period = 1800. ! en secondes
144	dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
145
146	OPEN (99, FILE='soil.def', STATUS='old', FORM='formatted', ERR=9999)
147	READ (99, *) min_period
148	READ (99, *) dalph_soil
149	PRINT *, 'Discretization for the soil model'
150	PRINT *, 'First level e-folding depth', min_period, ' dalph', &
151	dalph_soil
152	CLOSE (99)
153	9999 CONTINUE
154
155	! la premiere couche represente un dixieme de cycle diurne
156	fz1 = sqrt(min_period/3.14)
157
158	DO jk = 1, nsoilmx
159	rk1 = jk
160	rk2 = jk - 1
161	dz2(jk) = fz(rk1) - fz(rk2)
162	END DO
163	DO jk = 1, nsoilmx - 1
164	rk1 = jk + .5
165	rk2 = jk - .5
166	dz1(jk) = 1./(fz(rk1)-fz(rk2))
167	END DO
168	lambda = fz(.5)*dz1(1)
169	PRINT *, 'full layers, intermediate layers (seconds)'
170	DO jk = 1, nsoilmx
171	rk = jk
172	rk1 = jk + .5
173	rk2 = jk - .5
174	PRINT , 'fz=', fz(rk1)fz(rk2)3.14, fz(rk)fz(rk)*3.14
175	END DO
176	! PB
177	firstsurf(indice) = .FALSE.
178	! $$$ firstcall =.false.
179
180	! Initialisations:
181	! ----------------
182
183	ELSE !--not firstcall
184	! -----------------------------------------------------------------------
185	! Computation of the soil temperatures using the Cgrd and Dgrd
186	! coefficient computed at the previous time-step:
187	! -----------------------------------------------
188
189	! surface temperature
190	DO ig = 1, knon
191	ptsoil(ig, 1) = (lambdazc(ig,1,indice)+ptsrf(ig))/(lambda(1.-zd(ig,1, &
192	indice))+1.)
193	END DO
194
195	! other temperatures
196	DO jk = 1, nsoilmx - 1
197	DO ig = 1, knon
198	ptsoil(ig, jk+1) = zc(ig, jk, indice) + zd(ig, jk, indice)*ptsoil(ig, &
199	jk)
200	END DO
201	END DO
202
203	END IF !--not firstcall
204	! -----------------------------------------------------------------------
205	! Computation of the Cgrd and Dgrd coefficient for the next step:
206	! ---------------------------------------------------------------
207
208	! $$$ PB ajout pour cas glace de mer
209	IF (indice==is_sic) THEN
210	DO ig = 1, knon
211	ptsoil(ig, nsoilmx) = rtt - 1.8
212	END DO
213	END IF
214
215	DO jk = 1, nsoilmx
216	zdz2(jk) = dz2(jk)/ptimestep
217	END DO
218
219	DO ig = 1, knon
220	z1(ig, indice) = zdz2(nsoilmx) + dz1(nsoilmx-1)
221	zc(ig, nsoilmx-1, indice) = zdz2(nsoilmx)*ptsoil(ig, nsoilmx)/ &
222	z1(ig, indice)
223	zd(ig, nsoilmx-1, indice) = dz1(nsoilmx-1)/z1(ig, indice)
224	END DO
225
226	DO jk = nsoilmx - 1, 2, -1
227	DO ig = 1, knon
228	z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice)))
229	zc(ig, jk-1, indice) = (ptsoil(ig,jk)zdz2(jk)+dz1(jk)zc(ig,jk,indice) &
230	)*z1(ig, indice)
231	zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice)
232	END DO
233	END DO
234
235	! -----------------------------------------------------------------------
236	! computation of the surface diffusive flux from ground and
237	! calorific capacity of the ground:
238	! ---------------------------------
239
240	DO ig = 1, knon
241	pfluxgrd(ig) = ztherm_i(ig)dz1(1)(zc(ig,1,indice)+(zd(ig,1, &
242	indice)-1.)*ptsoil(ig,1))
243	pcapcal(ig) = ztherm_i(ig)(dz2(1)+ptimestep(1.-zd(ig,1,indice))*dz1(1))
244	z1(ig, indice) = lambda*(1.-zd(ig,1,indice)) + 1.
245	pcapcal(ig) = pcapcal(ig)/z1(ig, indice)
246	pfluxgrd(ig) = pfluxgrd(ig) + pcapcal(ig)(ptsoil(ig,1)z1(ig,indice)- &
247	lambda*zc(ig,1,indice)-ptsrf(ig))/ptimestep
248	END DO
249
250	contains
251
252	real function fz(rk)
253	real rk
254	fz = fz1(dalph_soil*rk-1.)/(dalph_soil-1.)
255	end function fz
256
257	END SUBROUTINE soil
258
259	end module soil_m