/[lmdze]/trunk/Sources/phylmd/Interface_surf/soil.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/Interface_surf/soil.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/soil.f revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC trunk/phylmd/soil.f revision 101 by guez, Mon Jul 7 17:45:21 2014 UTC
# Line 1  Line 1 
1  !  module soil_m
 ! $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,  
      s          pcapcal, pfluxgrd)  
       use dimens_m  
       use indicesol  
       use dimphy  
       use dimsoil  
       use SUPHEC_M  
       IMPLICIT NONE  
   
 c=======================================================================  
 c  
 c   Auteur:  Frederic Hourdin     30/01/92  
 c   -------  
 c  
 c   objet:  computation of : the soil temperature evolution  
 c   ------                   the surfacic heat capacity "Capcal"  
 c                            the surface conduction flux pcapcal  
 c  
 c  
 c   Method: implicit time integration  
 c   -------  
 c   Consecutive ground temperatures are related by:  
 c           T(k+1) = C(k) + D(k)*T(k)  (1)  
 c   the coefficients C and D are computed at the t-dt time-step.  
 c   Routine structure:  
 c   1)new temperatures are computed  using (1)  
 c   2)C and D coefficients are computed from the new temperature  
 c     profile for the t+dt time-step  
 c   3)the coefficients A and B are computed where the diffusive  
 c     fluxes at the t+dt time-step is given by  
 c            Fdiff = A + B Ts(t+dt)  
 c     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt  
 c            with F0 = A + B (Ts(t))  
 c                 Capcal = B*dt  
 c            
 c   Interface:  
 c   ----------  
 c  
 c   Arguments:  
 c   ----------  
 c   ptimestep            physical timestep (s)  
 c   indice               sub-surface index  
 c   snow(klon,nbsrf)     snow  
 c   ptsrf(klon)          surface temperature at time-step t (K)  
 c   ptsoil(klon,nsoilmx) temperature inside the ground (K)  
 c   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)  
 c   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)  
 c    
 c=======================================================================  
 c   declarations:  
 c   -------------  
   
   
 c-----------------------------------------------------------------------  
 c  arguments  
 c  ---------  
   
       REAL ptimestep  
       INTEGER indice, knon  
       REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon)  
       REAL pcapcal(klon),pfluxgrd(klon)  
   
 c-----------------------------------------------------------------------  
 c  local arrays  
 c  ------------  
   
       INTEGER ig,jk  
 c$$$      REAL zdz2(nsoilmx),z1(klon)  
       REAL zdz2(nsoilmx),z1(klon,nbsrf)  
       REAL min_period,dalph_soil  
       REAL ztherm_i(klon)  
   
 c   local saved variables:  
 c   ----------------------  
       REAL dz1(nsoilmx),dz2(nsoilmx)  
 c$$$          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./  
   
 c-----------------------------------------------------------------------  
 c   Depthts:  
 c   --------  
   
       REAL fz,rk,fz1,rk1,rk2  
       fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)  
       pfluxgrd(:) = 0.  
 c   calcul de l'inertie thermique a partir de la variable rnat.  
 c   on initialise a iice meme au-dessus d'un point de mer au cas  
 c   ou le point de mer devienne point de glace au pas suivant  
 c   on corrige si on a un point de terre avec ou sans glace  
 c  
       IF (indice.EQ.is_sic) THEN  
          DO ig = 1, knon  
             ztherm_i(ig)   = iice  
             IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno  
          ENDDO  
       ELSE IF (indice.EQ.is_lic) THEN  
          DO ig = 1, knon  
             ztherm_i(ig)   = iice  
             IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno  
          ENDDO  
       ELSE IF (indice.EQ.is_ter) THEN  
          DO ig = 1, knon  
             ztherm_i(ig)   = isol  
             IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno  
          ENDDO  
       ELSE IF (indice.EQ.is_oce) THEN  
          DO ig = 1, knon  
             ztherm_i(ig)   = iice  
          ENDDO  
       ELSE  
          PRINT*, "valeur d indice non prevue", indice  
          stop 1  
       ENDIF  
   
   
 c$$$      IF (firstcall) THEN  
       IF (firstsurf(indice)) THEN  
   
 c-----------------------------------------------------------------------  
 c   ground levels  
 c   grnd=z/l where l is the skin depth of the diurnal cycle:  
 c   --------------------------------------------------------  
   
          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,  
      s   '   dalph',dalph_soil  
          CLOSE(99)  
 9999     CONTINUE  
   
 c   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)  
          ENDDO  
          DO jk=1,nsoilmx-1  
             rk1=jk+.5  
             rk2=jk-.5  
             dz1(jk)=1./(fz(rk1)-fz(rk2))  
          ENDDO  
          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  
          ENDDO  
 C PB  
          firstsurf(indice) = .FALSE.  
 c$$$         firstcall =.false.  
   
 c   Initialisations:  
 c   ----------------  
   
       ELSE   !--not firstcall  
 c-----------------------------------------------------------------------  
 c   Computation of the soil temperatures using the Cgrd and Dgrd  
 c  coefficient computed at the previous time-step:  
 c  -----------------------------------------------  
   
 c    surface temperature  
          DO ig=1,knon  
             ptsoil(ig,1)=(lambda*zc(ig,1,indice)+ptsrf(ig))/  
      s      (lambda*(1.-zd(ig,1,indice))+1.)  
          ENDDO  
   
 c   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)  
             ENDDO  
          ENDDO  
   
       ENDIF !--not firstcall  
 c-----------------------------------------------------------------------  
 c   Computation of the Cgrd and Dgrd coefficient for the next step:  
 c   ---------------------------------------------------------------  
   
 c$$$  PB ajout pour cas glace de mer  
       IF (indice .EQ. is_sic) THEN  
           DO ig = 1 , knon  
             ptsoil(ig,nsoilmx) = RTT - 1.8  
           END DO  
       ENDIF  
   
       DO jk=1,nsoilmx  
          zdz2(jk)=dz2(jk)/ptimestep  
       ENDDO  
   
       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)  
       ENDDO  
   
       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)=  
      s      (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)  
          ENDDO  
       ENDDO  
   
 c-----------------------------------------------------------------------  
 c   computation of the surface diffusive flux from ground and  
 c   calorific capacity of the ground:  
 c   ---------------------------------  
   
       DO ig=1,knon  
          pfluxgrd(ig)=ztherm_i(ig)*dz1(1)*  
      s   (zc(ig,1,indice)+(zd(ig,1,indice)-1.)*ptsoil(ig,1))  
          pcapcal(ig)=ztherm_i(ig)*  
      s   (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)  
      s   + pcapcal(ig) * (ptsoil(ig,1) * z1(ig,indice)  
      $       - lambda * zc(ig,1,indice)  
      $       - ptsrf(ig))  
      s   /ptimestep  
       ENDDO  
2    
3        RETURN    IMPLICIT NONE
4        END  
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(klon)          surface temperature at time-step t (K)
53        ! ptsoil(klon,nsoilmx) temperature inside the ground (K)
54        ! pcapcal(klon)        surfacic specific heat (W*m-2*s*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(klon), 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 fz, rk, fz1, rk1, rk2
103    
104        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    
111        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    
135    
136        ! $$$      IF (firstcall) THEN
137        IF (firstsurf(indice)) THEN
138    
139           ! -----------------------------------------------------------------------
140           ! ground levels
141           ! grnd=z/l where l is the skin depth of the diurnal cycle:
142           ! --------------------------------------------------------
143    
144           min_period = 1800. ! en secondes
145           dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
146    
147           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) = (lambda*zc(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        ! -----------------------------------------------------------------------
206        ! Computation of the Cgrd and Dgrd coefficient for the next step:
207        ! ---------------------------------------------------------------
208    
209        ! $$$  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    
216        DO jk = 1, nsoilmx
217           zdz2(jk) = dz2(jk)/ptimestep
218        END DO
219    
220        DO ig = 1, knon
221           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        END DO
226    
227        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    
236        ! -----------------------------------------------------------------------
237        ! computation of the surface diffusive flux from ground and
238        ! calorific capacity of the ground:
239        ! ---------------------------------
240    
241        DO ig = 1, knon
242           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        END DO
250    
251      END SUBROUTINE soil
252    
253    end module soil_m

Legend:
Removed from v.38  
changed lines
  Added in v.101

  ViewVC Help
Powered by ViewVC 1.1.21