Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/hydrol.F90 @ 227

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

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

Download in other formats: