Context Navigation

hydrol.F90 @ 867

Last change on this file since 867 was 863, checked in by jleconte, 12 years ago

23/01/2013 == JL

Correction in largescale. a rneb factor was forgotten
Added some spectra in ave_stelpec
Corrected reevaporation in rain. Now conserve water better

Property svn:executable set to *

File size: 10.1 KB

Line
1	subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf, &
2	qsurf,dqsurf,dqs_hyd,pcapcal, &
3	albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
4
5	use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
6	USE surfdat_h
7	use comdiurn_h
8	USE comgeomfi_h
9	USE tracer_h
10
11	implicit none
12
13	!==================================================================
14	!
15	! Purpose
16	! -------
17	! Calculate the surface hydrology and albedo changes.
18	!
19	! Authors
20	! -------
21	! Adapted from LMDTERRE by B. Charnay (2010). Further
22	! modifications by R. Wordsworth (2010).
23	!
24	! Called by
25	! ---------
26	! physiq.F
27	!
28	! Calls
29	! -----
30	! none
31	!
32	! Notes
33	! -----
34	! rnat is terrain type: 0-ocean; 1-continent
35	!
36	!==================================================================
37
38	#include "dimensions.h"
39	#include "dimphys.h"
40	#include "comcstfi.h"
41	#include "callkeys.h"
42
43	integer ngrid,nq
44
45	! Inputs
46	! ------
47	real albedoocean
48	parameter (albedoocean=0.07)
49	real albedoice
50	save albedoice
51
52	real snowlayer
53	parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
54	real oceantime
55	parameter (oceantime=10243600)
56
57	logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
58	logical activerunoff ! enable simple runoff scheme?
59	logical oceanalbvary ! ocean albedo varies with the diurnal cycle?
60
61	! Arguments
62	! ---------
63	integer rnat(ngrid) ! I changed this to integer (RW)
64	real,dimension(:),allocatable,save :: runoff
65	real totalrunoff, tsea, oceanarea
66	save oceanarea
67
68	real ptimestep
69	real mu0(ngrid)
70	real qsurf(ngrid,nq), tsurf(ngrid)
71	real dqsurf(ngrid,nq), pdtsurf(ngrid)
72	real hice(ngrid)
73	real albedo0(ngrid), albedo(ngrid)
74
75	real oceanarea2
76
77	! Output
78	! ------
79	real dqs_hyd(ngrid,nq)
80	real pdtsurf_hyd(ngrid)
81
82	! Local
83	! -----
84	real a,b,E
85	integer ig,iq, icap ! wld like to remove icap
86	real fsnoi, subli, fauxo
87	real twater(ngrid)
88	real pcapcal(ngrid)
89	real hicebis(ngrid)
90	real zqsurf(ngrid,nq)
91	real ztsurf(ngrid)
92
93	integer, save :: ivap, iliq, iice
94
95	logical, save :: firstcall
96
97	data firstcall /.true./
98
99	oceanbulkavg=.false.
100	activerunoff=.false.
101	oceanalbvary=.false.
102
103	if(firstcall)then
104
105	ALLOCATE(runoff(ngrid))
106
107	ivap=igcm_h2o_vap
108	iliq=igcm_h2o_vap
109	iice=igcm_h2o_ice
110
111	write(,) "hydrol: ivap=",ivap
112	write(,) " iliq=",iliq
113	write(,) " iice=",iice
114
115	! Here's the deal: iice is used in place of igcm_h2o_ice both on the
116	! surface and in the atmosphere. ivap is used in
117	! place of igcm_h2o_vap ONLY in the atmosphere, while
118	! iliq is used in place of igcm_h2o_vap ONLY on the
119	! surface.
120
121	! Soon to be extended to the entire water cycle...
122
123	! Ice albedo = snow albedo for now
124	albedoice=albedosnow
125
126	! Total ocean surface area
127	oceanarea=0.
128	do ig=1,ngrid
129	if(rnat(ig).eq.0)then
130	oceanarea=oceanarea+area(ig)
131	endif
132	enddo
133
134	if(oceanbulkavg.and.(oceanarea.le.0.))then
135	print*,'How are we supposed to average the ocean'
136	print*,'temperature, when there are no oceans?'
137	call abort
138	endif
139
140	if(activerunoff.and.(oceanarea.le.0.))then
141	print*,'You have enabled runoff, but you have no oceans.'
142	print*,'Where did you think the water was going to go?'
143	call abort
144	endif
145
146	firstcall = .false.
147	endif
148
149	! add physical tendencies already calculated
150	! ------------------------------------------
151
152	do ig=1,ngrid
153	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
154	pdtsurf_hyd(ig)=0.0
155	do iq=1,nq
156	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
157	enddo
158	enddo
159
160	do ig=1,ngrid
161	do iq=1,nq
162	dqs_hyd(ig,iq) = 0.0
163	enddo
164	enddo
165
166	do ig = 1, ngrid
167
168	! Ocean
169	! -----
170	if(rnat(ig).eq.0)then
171
172	! re-calculate oceanic albedo
173	if(diurnal.and.oceanalbvary)then
174	fauxo = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
175	albedo(ig) = 1.1( .03 + .630/( 1. + fauxofauxo))
176	albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
177	else
178	albedo(ig) = albedoocean ! modif Benjamin
179	end if
180
181	! calculate oceanic ice height including the latent heat of ice formation
182	! hice is the height of oceanic ice with a maximum of maxicethick.
183	hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
184
185	! twater(ig) = tsurf(ig) + ptimestep*zdtsurf(ig) &
186	twater(ig) = ztsurf(ig) - hice(ig)RLFTTrhowater/pcapcal(ig)
187	! this is temperature water would have if we melted entire ocean ice layer
188	hicebis(ig) = hice(ig)
189	hice(ig) = 0.
190
191	if(twater(ig) .lt. T_h2O_ice_liq)then
192	E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))pcapcal(ig),RLFTTrhowater*maxicethick)
193	hice(ig) = E/(RLFTT*rhowater)
194	hice(ig) = max(hice(ig),0.0)
195	hice(ig) = min(hice(ig),maxicethick)
196	pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))RLFTTrhowater/pcapcal(ig)/ptimestep
197	albedo(ig) = albedoice
198
199	! if (zqsurf(ig,iice).ge.snowlayer) then
200	! albedo(ig) = albedoice
201	! else
202	! albedo(ig) = albedoocean &
203	! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
204	! endif
205
206	else
207
208	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
209	albedo(ig) = albedoocean
210
211	endif
212
213	zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
214	zqsurf(ig,iice) = hice(ig)*rhowater
215
216
217	! Continent
218	! ---------
219	elseif (rnat(ig).eq.1) then
220
221	! melt the snow
222	if(ztsurf(ig).gt.T_h2O_ice_liq)then
223	if(zqsurf(ig,iice).gt.1.0e-8)then
224
225	a = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
226	b = zqsurf(ig,iice)
227	fsnoi = min(a,b)
228
229	zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
230	zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi
231
232	! thermal effects
233	pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep
234
235	endif
236	else
237
238	! freeze the water
239	if(zqsurf(ig,iliq).gt.1.0e-8)then
240
241	a = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
242	b = zqsurf(ig,iliq)
243
244	fsnoi = min(a,b)
245
246	zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
247	zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi
248
249	! thermal effects
250	pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep
251
252	endif
253	endif
254
255	! deal with runoff
256	if(activerunoff)then
257
258	runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
259	if(ngrid.gt.1)then ! runoff only exists in 3D
260	if(runoff(ig).ne.0.0)then
261	zqsurf(ig,iliq) = mx_eau_sol
262	! runoff is added to ocean at end
263	endif
264	end if
265
266	endif
267
268	! re-calculate continental albedo
269	albedo(ig) = albedo0(ig) ! albedo0 = base values
270	if (zqsurf(ig,iice).ge.snowlayer) then
271	albedo(ig) = albedosnow
272	else
273	albedo(ig) = albedo0(ig) &
274	+ (albedosnow - albedo0(ig))*zqsurf(ig,iice)/snowlayer
275	endif
276
277	else
278
279	print*,'Surface type not recognised in hydrol.F!'
280	print*,'Exiting...'
281	call abort
282
283	endif
284
285	end do ! ig=1,ngrid
286
287	! perform crude bulk averaging of temperature in ocean
288	! ----------------------------------------------------
289	if(oceanbulkavg)then
290
291	oceanarea2=0.
292	DO ig=1,ngrid
293	if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
294	oceanarea2=oceanarea2+area(ig)*pcapcal(ig)
295	end if
296	END DO
297
298	tsea=0.
299	DO ig=1,ngrid
300	if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
301	tsea=tsea+ztsurf(ig)area(ig)pcapcal(ig)/oceanarea2
302	end if
303	END DO
304
305	DO ig=1,ngrid
306	if((rnat(ig).eq.0).and.(hice(ig).eq.0))then
307	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
308	end if
309	END DO
310
311	print*,'Mean ocean temperature = ',tsea,' K'
312
313	endif
314
315	! shove all the runoff water into the ocean
316	! -----------------------------------------
317	if(activerunoff)then
318
319	totalrunoff=0.
320	do ig=1,ngrid
321	if (rnat(ig).eq.1) then
322	totalrunoff = totalrunoff + area(ig)*runoff(ig)
323	endif
324	enddo
325
326	do ig=1,ngrid
327	if (rnat(ig).eq.0) then
328	zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
329	totalrunoff/oceanarea
330	endif
331	enddo
332
333	endif
334
335
336	! Re-add the albedo effects of CO2 ice if necessary
337	! -------------------------------------------------
338	if(co2cond)then
339
340	icap=1
341	do ig=1,ngrid
342	if (qsurf(ig,igcm_co2_ice).gt.0) then
343	albedo(ig) = albedice(icap)
344	endif
345	enddo
346
347	endif
348
349
350	do ig=1,ngrid
351	dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep
352	dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep
353	enddo
354
355	call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
356
357	return
358	end subroutine hydrol

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

Download in other formats:

Original Format