Context Navigation

soil.F @ 1526

Last change on this file since 1526 was 1524, checked in by emillour, 9 years ago

All GCMS:
More updates to enforce dynamics/physics separation:

get rid of references to "temps_mod" from physics packages;
make a "time_phylmdz_mod.F90" module to store that
information and fill it via "iniphysiq".

File size: 8.2 KB

Line
1	subroutine soil(ngrid,nsoil,firstcall,lastcall,
2	& therm_i,
3	& timestep,tsurf,tsoil,
4	& capcal,fluxgrd)
5
6	use comsoil_h, only: layer, mlayer, volcapa, inertiedat
7	use comcstfi_mod, only: pi
8	use time_phylmdz_mod, only: daysec
9	use planete_mod, only: year_day
10	use comgeomfi_h, only: long, lati
11
12	implicit none
13
14	!-----------------------------------------------------------------------
15	! Author: Ehouarn Millour
16	!
17	! Purpose: Compute soil temperature using an implict 1st order scheme
18	!
19	! Note: depths of layers and mid-layers, soil thermal inertia and
20	! heat capacity are commons in comsoil.h
21	!-----------------------------------------------------------------------
22
23	! include "dimensions.h"
24	! include "dimphys.h"
25
26
27	c-----------------------------------------------------------------------
28	! arguments
29	! ---------
30	! inputs:
31	integer,intent(in) :: ngrid ! number of (horizontal) grid-points
32	integer,intent(in) :: nsoil ! number of soil layers
33	logical,intent(in) :: firstcall ! identifier for initialization call
34	logical,intent(in) :: lastcall
35	real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia
36	real,intent(in) :: timestep ! time step
37	real,intent(in) :: tsurf(ngrid) ! surface temperature
38	! outputs:
39	real,intent(out) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
40	real,intent(out) :: capcal(ngrid) ! surface specific heat
41	real,intent(out) :: fluxgrd(ngrid) ! surface diffusive heat flux
42
43	! local saved variables:
44	real,dimension(:,:),save,allocatable :: mthermdiff ! mid-layer thermal diffusivity
45	real,dimension(:,:),save,allocatable :: thermdiff ! inter-layer thermal diffusivity
46	real,dimension(:),save,allocatable :: coefq ! q_{k+1/2} coefficients
47	real,dimension(:,:),save,allocatable :: coefd ! d_k coefficients
48	real,dimension(:,:),save,allocatable :: alph ! alpha_k coefficients
49	real,dimension(:,:),save,allocatable :: beta ! beta_k coefficients
50	real,save :: mu
51	!$OMP THREADPRIVATE(mthermdiff,thermdiff,coefq,coefd,alph,beta,mu)
52
53	! local variables:
54	integer ig,ik
55	real :: inertia_min,inertia_max
56	real :: diurnal_skin ! diurnal skin depth (m)
57	real :: annual_skin ! anuual skin depth (m)
58
59	! 0. Initialisations and preprocessing step
60	if (firstcall) then
61	! note: firstcall is set to .true. or .false. by the caller
62	! and not changed by soil.F
63
64	ALLOCATE(mthermdiff(ngrid,0:nsoil-1)) ! mid-layer thermal diffusivity
65	ALLOCATE(thermdiff(ngrid,nsoil-1)) ! inter-layer thermal diffusivity
66	ALLOCATE(coefq(0:nsoil-1)) ! q_{k+1/2} coefficients
67	ALLOCATE(coefd(ngrid,nsoil-1)) ! d_k coefficients
68	ALLOCATE(alph(ngrid,nsoil-1)) ! alpha_k coefficients
69	ALLOCATE(beta(ngrid,nsoil-1)) ! beta_k coefficients
70
71	! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
72	do ig=1,ngrid
73	do ik=0,nsoil-1
74	mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
75	! write(,),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik)
76	enddo
77	enddo
78
79	! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
80	do ig=1,ngrid
81	do ik=1,nsoil-1
82	thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)
83	& +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))
84	& /(mlayer(ik)-mlayer(ik-1))
85	! write(,),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
86	enddo
87	enddo
88
89	! 0.3 Build coefficients mu, q_{k+1/2}, d_k, alpha_k and capcal
90	! mu
91	mu=mlayer(0)/(mlayer(1)-mlayer(0))
92
93	! q_{1/2}
94	coefq(0)=volcapa*layer(1)/timestep
95	! q_{k+1/2}
96	do ik=1,nsoil-1
97	coefq(ik)=volcapa*(layer(ik+1)-layer(ik))
98	& /timestep
99	enddo
100
101	do ig=1,ngrid
102	! d_k
103	do ik=1,nsoil-1
104	coefd(ig,ik)=thermdiff(ig,ik)/(mlayer(ik)-mlayer(ik-1))
105	enddo
106
107	! alph_{N-1}
108	alph(ig,nsoil-1)=coefd(ig,nsoil-1)/
109	& (coefq(nsoil-1)+coefd(ig,nsoil-1))
110	! alph_k
111	do ik=nsoil-2,1,-1
112	alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*
113	& (1.-alph(ig,ik+1))+coefd(ig,ik))
114	enddo
115
116	! capcal
117	! Cstar
118	capcal(ig)=volcapa*layer(1)+
119	& (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
120	& (timestep*(1.-alph(ig,1)))
121	! Cs
122	capcal(ig)=capcal(ig)/(1.+mu(1.0-alph(ig,1))
123	& thermdiff(ig,1)/mthermdiff(ig,0))
124	!write(,)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
125	enddo ! of do ig=1,ngrid
126
127	! Additional checks: is the vertical discretization sufficient
128	! to resolve diurnal and annual waves?
129	do ig=1,ngrid
130	! extreme inertia for this column
131	inertia_min=minval(inertiedat(ig,:))
132	inertia_max=maxval(inertiedat(ig,:))
133	! diurnal and annual skin depth
134	diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
135	annual_skin=(inertia_max/volcapa)sqrt(year_daydaysec/pi)
136	if (0.5*diurnal_skin<layer(1)) then
137	! one should have the fist layer be at least half of diurnal skin depth
138	write(,) "soil Error: grid point ig=",ig
139	write(,) " longitude=",long(ig)*(180./pi)
140	write(,) " latitude=",lati(ig)*(180./pi)
141	write(,) " first soil layer depth ",layer(1)
142	write(,) " not small enough for a diurnal skin depth of ",
143	& diurnal_skin
144	write(,) " change soil layer distribution (comsoil_h.F90)"
145	stop
146	endif
147	if (2.*annual_skin>layer(nsoil)) then
148	! one should have the full soil be at least twice the diurnal skin depth
149	write(,) "soil Error: grid point ig=",ig
150	write(,) " longitude=",long(ig)*(180./pi)
151	write(,) " latitude=",lati(ig)*(180./pi)
152	write(,) " total soil layer depth ",layer(nsoil)
153	write(,) " not large enough for an annual skin depth of ",
154	& annual_skin
155	write(,) " change soil layer distribution (comsoil_h.F90)"
156	stop
157	endif
158	enddo ! of do ig=1,ngrid
159
160	else ! of if (firstcall)
161
162
163	! 1. Compute soil temperatures
164	! First layer:
165	do ig=1,ngrid
166	tsoil(ig,1)=(tsurf(ig)+mubeta(ig,1)
167	& thermdiff(ig,1)/mthermdiff(ig,0))/
168	& (1.+mu(1.0-alph(ig,1))
169	& thermdiff(ig,1)/mthermdiff(ig,0))
170	enddo
171	! Other layers:
172	do ik=1,nsoil-1
173	do ig=1,ngrid
174	tsoil(ig,ik+1)=alph(ig,ik)*tsoil(ig,ik)+beta(ig,ik)
175	enddo
176	enddo
177
178	endif! of if (firstcall)
179
180	! 2. Compute beta coefficients (preprocessing for next time step)
181	! Bottom layer, beta_{N-1}
182	do ig=1,ngrid
183	beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil)
184	& /(coefq(nsoil-1)+coefd(ig,nsoil-1))
185	enddo
186	! Other layers
187	do ik=nsoil-2,1,-1
188	do ig=1,ngrid
189	beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+
190	& coefd(ig,ik+1)*beta(ig,ik+1))/
191	& (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1))
192	& +coefd(ig,ik))
193	enddo
194	enddo
195
196
197	! 3. Compute surface diffusive flux & calorific capacity
198	do ig=1,ngrid
199	! Cstar
200	! capcal(ig)=volcapa(ig,1)*layer(ig,1)+
201	! & (thermdiff(ig,1)/(mlayer(ig,1)-mlayer(ig,0)))*
202	! & (timestep*(1.-alph(ig,1)))
203	! Fstar
204
205	! print*,'this far in soil 1'
206	! print*,'thermdiff=',thermdiff(ig,1)
207	! print*,'mlayer=',mlayer
208	! print*,'beta=',beta(ig,1)
209	! print*,'alph=',alph(ig,1)
210	! print*,'tsoil=',tsoil(ig,1)
211
212	fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
213	& (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))
214
215	! mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
216	! capcal(ig)=capcal(ig)/(1.+mu(1.0-alph(ig,1))
217	! & thermdiff(ig,1)/mthermdiff(ig,0))
218	! Fs
219	fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*
220	& (tsoil(ig,1)(1.+mu(1.0-alph(ig,1))*
221	& thermdiff(ig,1)/mthermdiff(ig,0))
222	& -tsurf(ig)-mubeta(ig,1)
223	& thermdiff(ig,1)/mthermdiff(ig,0))
224	enddo
225
226	end
227

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format