/[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/phylmd/soil.f revision 80 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/phylmd/soil.f90 revision 81 by guez, Wed Mar 5 14:38:41 2014 UTC
# Line 1  Line 1 
 !  
 ! $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  
1    
2        RETURN  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/soil.F,v 1.1.1.1 2004/05/19
3        END  ! 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 (W*m-2*s*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) = (lambda*zc(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

Legend:
Removed from v.80  
changed lines
  Added in v.81

  ViewVC Help
Powered by ViewVC 1.1.21