phylmd/Interface_surf/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	guez	3
2	guez	81	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/soil.F,v 1.1.1.1 2004/05/19
3			! 12:53:09 lmdzadmin Exp $
4	guez	3
5	guez	81	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	guez	3
14	guez	81	! =======================================================================
15	guez	3
16	guez	81	! Auteur: Frederic Hourdin 30/01/92
17			! -------
18	guez	3
19	guez	81	! objet: computation of : the soil temperature evolution
20			! ------ the surfacic heat capacity "Capcal"
21			! the surface conduction flux pcapcal
22	guez	3
23
24	guez	81	! 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	guez	3
40	guez	81	! Interface:
41			! ----------
42	guez	3
43	guez	81	! 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	guez	3
53	guez	81	! =======================================================================
54			! declarations:
55			! -------------
56	guez	3
57
58	guez	81	! -----------------------------------------------------------------------
59			! arguments
60			! ---------
61	guez	3
62	guez	81	REAL ptimestep
63			INTEGER indice, knon
64			REAL ptsrf(klon), ptsoil(klon, nsoilmx), snow(klon)
65			REAL pcapcal(klon), pfluxgrd(klon)
66	guez	3
67	guez	81	! -----------------------------------------------------------------------
68			! local arrays
69			! ------------
70	guez	3
71	guez	81	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	guez	3
77	guez	81	! 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	guez	3
89	guez	81	DATA firstcall/.TRUE./
90			DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./
91	guez	3
92	guez	81	DATA isol, isno, iice/2000., 2000., 2000./
93	guez	3
94	guez	81	! -----------------------------------------------------------------------
95			! Depthts:
96			! --------
97	guez	3
98	guez	81	REAL fz, rk, fz1, rk1, rk2
99	guez	3
100	guez	81	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	guez	3
107	guez	81	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	guez	3
131
132	guez	81	! $$$ IF (firstcall) THEN
133			IF (firstsurf(indice)) THEN
134	guez	3
135	guez	81	! -----------------------------------------------------------------------
136			! ground levels
137			! grnd=z/l where l is the skin depth of the diurnal cycle:
138			! --------------------------------------------------------
139	guez	3
140	guez	81	min_period = 1800. ! en secondes
141			dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
142	guez	3
143	guez	81	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	guez	3
152	guez	81	! la premiere couche represente un dixieme de cycle diurne
153			fz1 = sqrt(min_period/3.14)
154	guez	3
155	guez	81	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	guez	3
177	guez	81	! 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