Context Navigation

hydrol.F90 @ 1483

Last change on this file since 1483 was 1482, checked in by mturbet, 10 years ago
Implementation of the Spectral Albedo
Property svn:executable set to ``*
File size: 12.4 KB

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

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

Download in other formats:

Original Format