trunk/phylmd/soil.f


! $Header: /home/cvsroot/LMDZ4/libf/phylmd/soil.F,v 1.1.1.1 2004/05/19
! 12:53:09 lmdzadmin Exp $

SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil, pcapcal, &
    pfluxgrd)
  USE dimens_m
  USE indicesol
  USE dimphy
  USE dimsoil
  USE suphec_m
  IMPLICIT NONE

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

  ! 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

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