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" |
---|
23 | c----------------------------------------------------------------------- |
---|
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 | |
---|
134 | c ********* 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 |
---|
168 | c ****************************************** |
---|
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 | |
---|