Context Navigation

hydrol.F90 @ 1243

Last change on this file since 1243 was 965, checked in by emillour, 12 years ago

Common dynamics and generic/universal GCM:

LMDZ.COMMON: minor bug fix on the computation of physics mesh area in gcm.F
LMDZ.UNIVERSAL: missing clean initialization of tab_cntrl(:) array in phyredem.F90
LMDZ.GENERIC: minor bug fix in hydrol.F90, only output runoff if it is used. Update output routines so that all outputs files (stats, diagfi.nc, diagsoil.nc, diagspecIR.nc and diagspecVI.nc) can be generated when running LMDZ.UNIVERSAL in MPI mode.

Property svn:executable set to *

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

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

Download in other formats:

Original Format