Context Navigation

hydrol.F90 @ 747

Last change on this file since 747 was 650, checked in by jleconte, 13 years ago

Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
Minor name changes in watercommon
Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

Added consistency check in inifis
Moved 1d water initialization from physiqu to rcm1d

Property svn:executable set to *

File size: 10.2 KB

Line
1	subroutine hydrol(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
7	implicit none
8
9	!==================================================================
10	!
11	! Purpose
12	! -------
13	! Calculate the surface hydrology and albedo changes.
14	!
15	! Authors
16	! -------
17	! Adapted from LMDTERRE by B. Charnay (2010). Further
18	! modifications by R. Wordsworth (2010).
19	!
20	! Called by
21	! ---------
22	! physiq.F
23	!
24	! Calls
25	! -----
26	! none
27	!
28	! Notes
29	! -----
30	! rnat is terrain type: 0-ocean; 1-continent
31	!
32	!==================================================================
33
34	#include "dimensions.h"
35	#include "dimphys.h"
36	#include "comcstfi.h"
37	#include "callkeys.h"
38	#include "tracer.h"
39	#include "fisice.h"
40	#include "comgeomfi.h"
41	#include "comdiurn.h"
42	#include "surfdat.h"
43
44	! Inputs
45	! ------
46	real albedoocean
47	parameter (albedoocean=0.07)
48	real albedoice
49	save albedoice
50
51	real snowlayer
52	parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
53	real oceantime
54	parameter (oceantime=10243600)
55
56	logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
57	logical activerunoff ! enable simple runoff scheme?
58	logical oceanalbvary ! ocean albedo varies with the diurnal cycle?
59
60	! Arguments
61	! ---------
62	integer rnat(ngridmx) ! I changed this to integer (RW)
63	real runoff(ngridmx)
64	save runoff
65	real totalrunoff, tsea, oceanarea
66	save oceanarea
67
68	real ptimestep
69	real mu0(ngridmx)
70	real qsurf(ngridmx,nqmx), tsurf(ngridmx)
71	real dqsurf(ngridmx,nqmx), pdtsurf(ngridmx)
72	real hice(ngridmx)
73	real albedo0(ngridmx), albedo(ngridmx)
74
75	real oceanarea2
76
77	! Output
78	! ------
79	real dqs_hyd(ngridmx,nqmx)
80	real pdtsurf_hyd(ngridmx)
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(ngridmx)
88	real pcapcal(ngridmx)
89	real hicebis(ngridmx)
90	real zqsurf(ngridmx,nqmx)
91	real ztsurf(ngridmx)
92
93	integer ivap, iliq, iice
94	save ivap, iliq, iice
95
96	logical firstcall
97	save firstcall
98
99	data firstcall /.true./
100
101	oceanbulkavg=.false.
102	activerunoff=.false.
103	oceanalbvary=.false.
104
105	if(firstcall)then
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,ngridmx
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,ngridmx
153	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
154	pdtsurf_hyd(ig)=0.0
155	do iq=1,nqmx
156	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
157	enddo
158	enddo
159
160	do ig=1,ngridmx
161	do iq=1,nqmx
162	dqs_hyd(ig,iq) = 0.0
163	enddo
164	enddo
165
166	do ig = 1, ngridmx
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))RLFTT &
197	rhowater/pcapcal(ig)/ptimestep
198	albedo(ig) = albedoice
199
200	! if (zqsurf(ig,iice).ge.snowlayer) then
201	! albedo(ig) = albedoice
202	! else
203	! albedo(ig) = albedoocean &
204	! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
205	! endif
206
207	else
208
209	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
210	albedo(ig) = albedoocean
211
212	endif
213
214	zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
215	zqsurf(ig,iice) = hice(ig)*rhowater
216
217
218	! Continent
219	! ---------
220	elseif (rnat(ig).eq.1) then
221
222	! melt the snow
223	if(ztsurf(ig).gt.T_h2O_ice_liq)then
224	if(zqsurf(ig,iice).gt.1.0e-8)then
225
226	a = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
227	b = zqsurf(ig,iice)
228	fsnoi = min(a,b)
229
230	zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
231	zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi
232
233	! thermal effects
234	pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep
235
236	endif
237	else
238
239	! freeze the water
240	if(zqsurf(ig,iliq).gt.1.0e-8)then
241
242	a = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
243	b = zqsurf(ig,iliq)
244
245	fsnoi = min(a,b)
246
247	zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
248	zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi
249
250	! thermal effects
251	pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep
252
253	endif
254	endif
255
256	! deal with runoff
257	if(activerunoff)then
258
259	runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
260	if(ngridmx.gt.1)then ! runoff only exists in 3D
261	if(runoff(ig).ne.0.0)then
262	zqsurf(ig,iliq) = mx_eau_sol
263	! runoff is added to ocean at end
264	endif
265	end if
266
267	endif
268
269	! re-calculate continental albedo
270	albedo(ig) = albedo0(ig) ! albedo0 = base values
271	if (zqsurf(ig,iice).ge.snowlayer) then
272	albedo(ig) = albedosnow
273	else
274	albedo(ig) = albedo0(ig) &
275	+ (albedosnow - albedo0(ig))*zqsurf(ig,iice)/snowlayer
276	endif
277
278	else
279
280	print*,'Surface type not recognised in hydrol.F!'
281	print*,'Exiting...'
282	call abort
283
284	endif
285
286	end do ! ig=1,ngridmx
287
288	! perform crude bulk averaging of temperature in ocean
289	! ----------------------------------------------------
290	if(oceanbulkavg)then
291
292	oceanarea2=0.
293	DO ig=1,ngridmx
294	if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
295	oceanarea2=oceanarea2+area(ig)*pcapcal(ig)
296	end if
297	END DO
298
299	tsea=0.
300	DO ig=1,ngridmx
301	if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
302	tsea=tsea+ztsurf(ig)area(ig)pcapcal(ig)/oceanarea2
303	end if
304	END DO
305
306	DO ig=1,ngridmx
307	if((rnat(ig).eq.0).and.(hice(ig).eq.0))then
308	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
309	end if
310	END DO
311
312	print*,'Mean ocean temperature = ',tsea,' K'
313
314	endif
315
316	! shove all the runoff water into the ocean
317	! -----------------------------------------
318	if(activerunoff)then
319
320	totalrunoff=0.
321	do ig=1,ngridmx
322	if (rnat(ig).eq.1) then
323	totalrunoff = totalrunoff + area(ig)*runoff(ig)
324	endif
325	enddo
326
327	do ig=1,ngridmx
328	if (rnat(ig).eq.0) then
329	zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
330	totalrunoff/oceanarea
331	endif
332	enddo
333
334	endif
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,ngridmx
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,ngridmx
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(ngridmx,'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