Context Navigation

soil.F @ 3523

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

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

Download in other formats:

Original Format