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	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			! 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	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			REAL ptsrf(klon), ptsoil(klon, nsoilmx), snow(klon)
69			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	101	REAL fz, rk, fz1, rk1, rk2
103	guez	3
104	guez	101	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	guez	3
111	guez	101	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	guez	3
135
136	guez	101	! $$$ IF (firstcall) THEN
137			IF (firstsurf(indice)) THEN
138	guez	3
139	guez	101	! -----------------------------------------------------------------------
140			! ground levels
141			! grnd=z/l where l is the skin depth of the diurnal cycle:
142			! --------------------------------------------------------
143	guez	3
144	guez	101	min_period = 1800. ! en secondes
145			dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
146	guez	81
147	guez	101	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	guez	81	! -----------------------------------------------------------------------
206	guez	101	! Computation of the Cgrd and Dgrd coefficient for the next step:
207			! ---------------------------------------------------------------
208	guez	81
209	guez	101	! $$$ 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	guez	81
216	guez	101	DO jk = 1, nsoilmx
217			zdz2(jk) = dz2(jk)/ptimestep
218	guez	81	END DO
219
220			DO ig = 1, knon
221	guez	101	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	guez	81	END DO
226
227	guez	101	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	guez	81
236	guez	101	! -----------------------------------------------------------------------
237			! computation of the surface diffusive flux from ground and
238			! calorific capacity of the ground:
239			! ---------------------------------
240	guez	81
241			DO ig = 1, knon
242	guez	101	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	guez	81	END DO
250
251	guez	101	END SUBROUTINE soil
252	guez	81
253	guez	101	end module soil_m