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	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	101	SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil, pcapcal, &
8			pfluxgrd)
9	guez	3
10	guez	101	! From LMDZ4/libf/phylmd/soil.F, version 1.1.1.1 2004/05/19
11	guez	3
12	guez	101	USE dimens_m
13			USE indicesol
14			USE dimphy
15			USE dimsoil
16			USE suphec_m
17	guez	3
18	guez	101	! =======================================================================
19	guez	3
20	guez	101	! Auteur: Frederic Hourdin 30/01/92
21			! -------
22	guez	3
23	guez	101	! objet: computation of : the soil temperature evolution
24			! ------ the surfacic heat capacity "Capcal"
25			! the surface conduction flux pcapcal
26	guez	3
27
28	guez	101	! 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	guez	3
44	guez	101	! Interface:
45			! ----------
46	guez	3
47	guez	101	! Arguments:
48			! ----------
49			! ptimestep physical timestep (s)
50			! indice sub-surface index
51			! snow(klon,nbsrf) snow
52	guez	104	! ptsrf(knon) surface temperature at time-step t (K)
53	guez	101	! 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	guez	3
57	guez	101	! =======================================================================
58			! declarations:
59			! -------------
60	guez	3
61
62	guez	101	! -----------------------------------------------------------------------
63			! arguments
64			! ---------
65	guez	3
66	guez	101	REAL ptimestep
67			INTEGER indice, knon
68	guez	104	REAL ptsrf(knon), ptsoil(klon, nsoilmx), snow(klon)
69	guez	101	REAL pcapcal(klon), pfluxgrd(klon)
70	guez	3
71	guez	101	! -----------------------------------------------------------------------
72			! local arrays
73			! ------------
74	guez	3
75	guez	101	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	guez	3
81	guez	101	! 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	guez	3
93	guez	101	DATA firstcall/.TRUE./
94			DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./
95	guez	3
96	guez	101	DATA isol, isno, iice/2000., 2000., 2000./
97	guez	3
98	guez	101	! -----------------------------------------------------------------------
99			! Depthts:
100			! --------
101	guez	3
102	guez	108	REAL rk, fz1, rk1, rk2
103	guez	3
104	guez	101	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	guez	3
110	guez	101	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	guez	3
134
135	guez	101	! $$$ IF (firstcall) THEN
136			IF (firstsurf(indice)) THEN
137	guez	3
138	guez	101	! -----------------------------------------------------------------------
139			! ground levels
140			! grnd=z/l where l is the skin depth of the diurnal cycle:
141			! --------------------------------------------------------
142	guez	3
143	guez	101	min_period = 1800. ! en secondes
144			dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
145	guez	81
146	guez	101	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	guez	81	! -----------------------------------------------------------------------
205	guez	101	! Computation of the Cgrd and Dgrd coefficient for the next step:
206			! ---------------------------------------------------------------
207	guez	81
208	guez	101	! $$$ 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	guez	81
215	guez	101	DO jk = 1, nsoilmx
216			zdz2(jk) = dz2(jk)/ptimestep
217	guez	81	END DO
218
219			DO ig = 1, knon
220	guez	101	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	guez	81	END DO
225
226	guez	101	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	guez	81
235	guez	101	! -----------------------------------------------------------------------
236			! computation of the surface diffusive flux from ground and
237			! calorific capacity of the ground:
238			! ---------------------------------
239	guez	81
240			DO ig = 1, knon
241	guez	101	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	guez	81	END DO
249
250	guez	108	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	guez	101	END SUBROUTINE soil
258	guez	81
259	guez	101	end module soil_m