source: trunk/LMDZ.PLUTO.old/libf/phypluto/soil.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM

File size: 7.2 KB
Line 
1      subroutine soil(ngrid,nsoil,firstcall,
2     &          therm_i,
3     &          timestep,tsurf,tsoil,
4     &          capcal,fluxgrd)
5
6      use comgeomfi_h
7      use comsoil_h
8      implicit none
9
10!-----------------------------------------------------------------------
11!  Author: Ehouarn Millour
12!
13!  Purpose: Compute soil temperature using an implict 1st order scheme
14
15!  Note: depths of layers and mid-layers, soil thermal inertia and
16!        heat capacity are commons in comsoil.h
17!-----------------------------------------------------------------------
18
19#include "dimensions.h"
20#include "dimphys.h"
21#include "callkeys.h"
22#include "comcstfi.h"
23c-----------------------------------------------------------------------
24!  arguments
25!  ---------
26!  inputs:
27      integer ngrid     ! number of (horizontal) grid-points
28      integer nsoil     ! number of soil layers
29      logical firstcall ! identifier for initialization call
30      real therm_i(ngrid,nsoil) ! thermal inertia
31      real timestep         ! time step
32      real tsurf(ngrid)   ! surface temperature
33! outputs:
34      real tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
35      real capcal(ngrid) ! surface specific heat
36      real fluxgrd(ngrid) ! surface diffusive heat flux
37
38! local saved variables:
39!      real,save :: layer(ngridmx,nsoilmx)      ! layer depth
40      real,save :: mthermdiff(ngridmx,0:nsoilmx-1) ! mid-layer thermal diffusivity
41      real,save :: thermdiff(ngridmx,nsoilmx-1) ! inter-layer thermal diffusivity
42      real,save :: coefq(0:nsoilmx-1)           ! q_{k+1/2} coefficients
43      real,save :: coefd(ngridmx,nsoilmx-1)     ! d_k coefficients
44      real,save :: alph(ngridmx,nsoilmx-1)      ! alpha_k coefficients
45      real,save :: beta(ngridmx,nsoilmx-1)      ! beta_k coefficients
46      real,save :: mu
47     
48! local variables:
49      integer ig,ik
50
51! 0. Initialisations and preprocessing step
52      if (firstcall.or.changeti) then
53      ! note: firstcall is set to .true. or .false. by the caller
54      !       and not changed by soil.F
55! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
56      do ig=1,ngrid
57        do ik=0,nsoil-1
58          !write(*,*),'soil: ik: ',ik,' TI',therm_i(ig,ik+1)
59          mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
60!         write(*,*),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik)
61        enddo
62      enddo
63
64! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
65      do ig=1,ngrid
66        do ik=1,nsoil-1
67      thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)
68     &                +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))
69     &                    /(mlayer(ik)-mlayer(ik-1))
70!       write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
71        enddo
72      enddo
73 
74! 0.3 Build coefficients mu, q_{k+1/2}, d_k, alpha_k and capcal
75      ! mu
76      mu=mlayer(0)/(mlayer(1)-mlayer(0))
77
78      ! q_{1/2}
79      coefq(0)=volcapa*layer(1)/timestep
80        ! q_{k+1/2}
81        do ik=1,nsoil-1
82          coefq(ik)=volcapa*(layer(ik+1)-layer(ik))
83     &                 /timestep
84        enddo
85
86      do ig=1,ngrid
87        ! d_k
88        do ik=1,nsoil-1
89          coefd(ig,ik)=thermdiff(ig,ik)/(mlayer(ik)-mlayer(ik-1))
90        enddo
91       
92        ! alph_{N-1}
93        alph(ig,nsoil-1)=coefd(ig,nsoil-1)/
94     &                  (coefq(nsoil-1)+coefd(ig,nsoil-1))
95        ! alph_k
96        do ik=nsoil-2,1,-1
97          alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*
98     &                              (1.-alph(ig,ik+1))+coefd(ig,ik))
99        enddo
100
101        ! capcal
102! Cstar
103        capcal(ig)=volcapa*layer(1)+
104     &              (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
105     &              (timestep*(1.-alph(ig,1)))
106! Cs
107        !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
108        !write(*,*)'thermdiff=',thermdiff
109        !write(*,*)'mthermdiff=',mthermdiff
110        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
111     &                         thermdiff(ig,1)/mthermdiff(ig,0))
112        !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
113      enddo ! of do ig=1,ngrid
114           
115      endif
116
117      if (.not.firstcall) then
118
119!  1. Compute soil temperatures
120! First layer:
121      do ig=1,ngrid
122        tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)*
123     &                         thermdiff(ig,1)/mthermdiff(ig,0))/
124     &              (1.+mu*(1.0-alph(ig,1))*
125     &               thermdiff(ig,1)/mthermdiff(ig,0))
126      enddo
127! Other layers:
128      do ik=1,nsoil-1
129        do ig=1,ngrid
130          tsoil(ig,ik+1)=alph(ig,ik)*tsoil(ig,ik)+beta(ig,ik)
131        enddo
132      enddo
133
134c  ********* Flux at the bottom
135      do ig=1,ngrid
136
137       if (assymflux) then
138         if ( (lati(ig)*180./pi.le.90.).and.
139     &           (lati(ig)*180./pi.ge.0.)  ) then
140             tsoil(ig,nsoil) = tsoil(ig,nsoil)
141     &       + timestep*fluxgeo2/(volcapa*(layer(nsoil)-layer(nsoil-1)))
142         else   
143             tsoil(ig,nsoil) = tsoil(ig,nsoil)
144     &       + timestep*fluxgeo/(volcapa*(layer(nsoil)-layer(nsoil-1)))
145         endif
146
147       elseif(patchflux.eq.1) then
148         if ( (long(ig)*180./pi.le.180.).and.(long(ig)*180./pi.ge.174.)
149     &   .and.(((lati(ig)*180./pi.le.46.).and.(lati(ig)*180./pi.ge.42.))
150     &   .or.  ((lati(ig)*180./pi.le.36.).and.(lati(ig)*180./pi.ge.32.))
151     &   .or.  ((lati(ig)*180./pi.le.26.).and.(lati(ig)*180./pi.ge.22.))
152     &   .or.  ((lati(ig)*180./pi.le.16.).and.(lati(ig)*180./pi.ge.12.))
153     &   .or.  ((lati(ig)*180./pi.le.6.).and.(lati(ig)*180./pi.ge.2.))
154     &           ) ) then
155             tsoil(ig,nsoil) = tsoil(ig,nsoil)
156     &       + timestep*fluxgeo2/(volcapa*(layer(nsoil)-layer(nsoil-1)))
157         else   
158             tsoil(ig,nsoil) = tsoil(ig,nsoil)
159     &       + timestep*fluxgeo/(volcapa*(layer(nsoil)-layer(nsoil-1)))
160         endif
161
162       else
163           tsoil(ig,nsoil) = tsoil(ig,nsoil)
164     &     + timestep*fluxgeo/(volcapa*(layer(nsoil)-layer(nsoil-1)))
165
166       endif
167      enddo
168c ******************************************
169     
170      endif! of if (not firstcall)
171
172!  2. Compute beta coefficients (preprocessing for next time step)
173! Bottom layer, beta_{N-1}
174      do ig=1,ngrid
175        beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil)
176     &                   /(coefq(nsoil-1)+coefd(ig,nsoil-1))
177      enddo
178! Other layers
179      do ik=nsoil-2,1,-1
180        do ig=1,ngrid
181          beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+
182     &                 coefd(ig,ik+1)*beta(ig,ik+1))/
183     &                 (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1))
184     &                  +coefd(ig,ik))
185        enddo
186      enddo
187
188!  3. Compute surface diffusive flux & calorific capacity
189      do ig=1,ngrid
190! Cstar
191!        capcal(ig)=volcapa(ig,1)*layer(ig,1)+
192!     &              (thermdiff(ig,1)/(mlayer(ig,1)-mlayer(ig,0)))*
193!     &              (timestep*(1.-alph(ig,1)))
194! Fstar
195        fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
196     &              (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))
197
198!        mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
199!        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
200!     &                         thermdiff(ig,1)/mthermdiff(ig,0))
201! Fs
202        fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*
203     &              (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))*
204     &                         thermdiff(ig,1)/mthermdiff(ig,0))
205     &               -tsurf(ig)-mu*beta(ig,1)*
206     &                          thermdiff(ig,1)/mthermdiff(ig,0))
207
208      enddo
209      end
210
Note: See TracBrowser for help on using the repository browser.