Context Navigation

hydrol.F90 @ 1477

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

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

Download in other formats:

Original Format