/[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/Interface_surf/soil.f revision 108 by guez, Tue Sep 16 14:00:41 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(knon)          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(knon), 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 rk, fz1, rk1, rk2
103    
104        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    
110        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    
134    
135        ! $$$      IF (firstcall) THEN
136        IF (firstsurf(indice)) THEN
137    
138           ! -----------------------------------------------------------------------
139           ! ground levels
140           ! grnd=z/l where l is the skin depth of the diurnal cycle:
141           ! --------------------------------------------------------
142    
143           min_period = 1800. ! en secondes
144           dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
145    
146           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) = (lambda*zc(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        ! -----------------------------------------------------------------------
205        ! Computation of the Cgrd and Dgrd coefficient for the next step:
206        ! ---------------------------------------------------------------
207    
208        ! $$$  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    
215        DO jk = 1, nsoilmx
216           zdz2(jk) = dz2(jk)/ptimestep
217        END DO
218    
219        DO ig = 1, knon
220           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        END DO
225    
226        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    
235        ! -----------------------------------------------------------------------
236        ! computation of the surface diffusive flux from ground and
237        ! calorific capacity of the ground:
238        ! ---------------------------------
239    
240        DO ig = 1, knon
241           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        END DO
249    
250      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      END SUBROUTINE soil
258    
259    end module soil_m
Legend:
Removed from v.38  
changed lines
  Added in v.108
  ViewVC Help
Powered by ViewVC 1.1.21