Context Navigation

hydrol.F90 @ 1521

Last change on this file since 1521 was 1521, checked in by emillour, 9 years ago

All GCMs: Updates to make planetary codes (+Earth) setups converge.

Made a "phy_common" directory in libf, to contain routines common (wrt structural nature of underlying code/grid) to all LMDZ-related physics packages.
moved all "mod_phys_*" and "mod_grid_phy_lmdz" files from dynlonlat_phylonlat to "phy_common"
moved "ioipsl_getincom_p.F90 from "misc" to "phy_common" and modified it to match Earth GCM version and renamed it ioipsl_getin_p_mod.F90
added an "abort_physics" (as in Earth GCM) in "phy_common"
added a "print_control_mod.F90 (as in Earth GCM) in phy_common
made similar changes in LMDZ.GENERIC and LMDZ.MARS

Property svn:executable set to *

File size: 12.5 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_getin_p_mod, only: getin_p
9	use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
10	USE surfdat_h
11	use comdiurn_h
12	USE comgeomfi_h
13	USE tracer_h
14	use slab_ice_h
15	use callkeys_mod, only: albedosnow,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
16	use radinc_h, only : L_NSPECTV
17
18	implicit none
19
20	!==================================================================
21	!
22	! Purpose
23	! -------
24	! Calculate the surface hydrology and albedo changes.
25	!
26	! Authors
27	! -------
28	! Adapted from LMDTERRE by B. Charnay (2010). Further
29	! Modifications by R. Wordsworth (2010).
30	! Spectral albedo by M. Turbet (2015).
31	!
32	! Called by
33	! ---------
34	! physiq.F
35	!
36	! Calls
37	! -----
38	! none
39	!
40	! Notes
41	! -----
42	! rnat is terrain type: 0-ocean; 1-continent
43	!
44	!==================================================================
45
46	integer ngrid,nq
47
48	! Inputs
49	! ------
50	real snowlayer
51	parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
52	real oceantime
53	parameter (oceantime=10243600)
54
55	logical,save :: oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
56	logical,save :: activerunoff ! enable simple runoff scheme?
57	logical,save :: oceanalbvary ! ocean albedo varies with the diurnal cycle?
58	!$OMP THREADPRIVATE(oceanbulkavg,activerunoff,oceanalbvary)
59
60	! Arguments
61	! ---------
62	real rnat(ngrid) ! I changed this to integer (RW)
63	real,dimension(:),allocatable,save :: runoff
64	real totalrunoff, tsea, oceanarea
65	save oceanarea
66	!$OMP THREADPRIVATE(runoff,oceanarea)
67
68	real ptimestep
69	real mu0(ngrid)
70	real qsurf(ngrid,nq), tsurf(ngrid)
71	real dqsurf(ngrid,nq), pdtsurf(ngrid)
72	real hice(ngrid)
73	real albedo(ngrid,L_NSPECTV)
74	real albedo_bareground(ngrid)
75	real albedo_snow_SPECTV(L_NSPECTV)
76	real albedo_co2_ice_SPECTV(L_NSPECTV)
77	real pctsrf_sic(ngrid), sea_ice(ngrid)
78
79	real oceanarea2
80
81	! Output
82	! ------
83	real dqs_hyd(ngrid,nq)
84	real pdtsurf_hyd(ngrid)
85
86	! Local
87	! -----
88	real a,b,E
89	integer ig,iq, nw
90	real fsnoi, subli, fauxo
91	real twater(ngrid)
92	real pcapcal(ngrid)
93	real hicebis(ngrid)
94	real zqsurf(ngrid,nq)
95	real ztsurf(ngrid)
96	real albedo_sic, alb_ice
97	real zfra
98
99	integer, save :: ivap, iliq, iice
100	!$OMP THREADPRIVATE(ivap,iliq,iice)
101
102	logical, save :: firstcall
103	!$OMP THREADPRIVATE(firstcall)
104
105	data firstcall /.true./
106
107
108	if(firstcall)then
109
110	oceanbulkavg=.false.
111	oceanalbvary=.false.
112	write(,)"Activate runnoff into oceans?"
113	activerunoff=.false.
114	call getin_p("activerunoff",activerunoff)
115	write(,)" activerunoff = ",activerunoff
116
117
118
119	if (activerunoff) ALLOCATE(runoff(ngrid))
120
121	ivap=igcm_h2o_vap
122	iliq=igcm_h2o_vap
123	iice=igcm_h2o_ice
124
125	write(,) "hydrol: ivap=",ivap
126	write(,) " iliq=",iliq
127	write(,) " iice=",iice
128
129	! Here's the deal: iice is used in place of igcm_h2o_ice both on the
130	! surface and in the atmosphere. ivap is used in
131	! place of igcm_h2o_vap ONLY in the atmosphere, while
132	! iliq is used in place of igcm_h2o_vap ONLY on the
133	! surface.
134	! Soon to be extended to the entire water cycle...
135
136	! Total ocean surface area
137	oceanarea=0.
138	do ig=1,ngrid
139	if(nint(rnat(ig)).eq.0)then
140	oceanarea=oceanarea+area(ig)
141	endif
142	enddo
143
144	if(oceanbulkavg.and.(oceanarea.le.0.))then
145	print*,'How are we supposed to average the ocean'
146	print*,'temperature, when there are no oceans?'
147	call abort
148	endif
149
150	if(activerunoff.and.(oceanarea.le.0.))then
151	print*,'You have enabled runoff, but you have no oceans.'
152	print*,'Where did you think the water was going to go?'
153	call abort
154	endif
155
156	firstcall = .false.
157	endif
158
159	! add physical tendencies already calculated
160	! ------------------------------------------
161
162	do ig=1,ngrid
163	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
164	pdtsurf_hyd(ig)=0.0
165	do iq=1,nq
166	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
167	enddo
168	enddo
169
170	do ig=1,ngrid
171	do iq=1,nq
172	dqs_hyd(ig,iq) = 0.0
173	enddo
174	enddo
175
176	do ig = 1, ngrid
177
178	! Ocean
179	! -----
180	if(nint(rnat(ig)).eq.0)then
181
182	! re-calculate oceanic albedo
183	! if(diurnal.and.oceanalbvary)then
184	! fauxo = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
185	! albedo(ig) = 1.1( .03 + .630/( 1. + fauxofauxo))
186	! albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
187	! else
188	do nw=1,L_NSPECTV
189	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.
190	enddo
191	! end if
192
193
194	if(ok_slab_ocean) then
195
196	zfra = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/45.0)) ! Snow Fraction (Critical height 45kg/m2~15cm)
197	alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) & ! Ice Albedo
198	*exp(-sea_ice(ig)/h_alb_ice)
199	! Albedo final calculation :
200	do nw=1,L_NSPECTV
201	albedo(ig,nw) = pctsrf_sic(ig)* &
202	(albedo_snow_SPECTV(nw)zfra + alb_ice(1.0-zfra)) &
203	+ (1.-pctsrf_sic(ig))*alb_ocean
204	enddo
205
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	do nw=1,L_NSPECTV
228	albedo(ig,nw) = albedo_snow_SPECTV(nw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
229	enddo
230
231	! if (zqsurf(ig,iice).ge.snowlayer) then
232	! albedo(ig) = albedoice
233	! else
234	! albedo(ig) = albedoocean &
235	! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
236	! endif
237
238	else
239
240	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
241	DO nw=1,L_NSPECTV
242	albedo(ig,nw) = alb_ocean
243	ENDDO
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	DO nw=1,L_NSPECTV
306	albedo(ig,nw) = albedo_bareground(ig)
307	ENDDO
308	if (zqsurf(ig,iice).ge.snowlayer) then
309	DO nw=1,L_NSPECTV
310	albedo(ig,nw) = albedo_snow_SPECTV(nw)
311	ENDDO
312	else
313	DO nw=1,L_NSPECTV
314	albedo(ig,nw) = albedo_bareground(ig) &
315	+ (albedo_snow_SPECTV(nw) - albedo_bareground(ig)) &
316	*zqsurf(ig,iice)/snowlayer
317	ENDDO
318	endif
319
320	else
321
322	print*,'Surface type not recognised in hydrol.F!'
323	print*,'Exiting...'
324	call abort
325
326	endif
327
328	end do ! ig=1,ngrid
329
330	! perform crude bulk averaging of temperature in ocean
331	! ----------------------------------------------------
332	if(oceanbulkavg)then
333
334	oceanarea2=0.
335	DO ig=1,ngrid
336	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
337	oceanarea2=oceanarea2+area(ig)*pcapcal(ig)
338	end if
339	END DO
340
341	tsea=0.
342	DO ig=1,ngrid
343	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
344	tsea=tsea+ztsurf(ig)area(ig)pcapcal(ig)/oceanarea2
345	end if
346	END DO
347
348	DO ig=1,ngrid
349	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0))then
350	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
351	end if
352	END DO
353
354	print*,'Mean ocean temperature = ',tsea,' K'
355
356	endif
357
358	! shove all the runoff water into the ocean
359	! -----------------------------------------
360	if(activerunoff)then
361
362	totalrunoff=0.
363	do ig=1,ngrid
364	if (nint(rnat(ig)).eq.1) then
365	totalrunoff = totalrunoff + area(ig)*runoff(ig)
366	endif
367	enddo
368
369	do ig=1,ngrid
370	if (nint(rnat(ig)).eq.0) then
371	zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
372	totalrunoff/oceanarea
373	endif
374	enddo
375
376	endif
377
378
379	! Re-add the albedo effects of CO2 ice if necessary
380	! -------------------------------------------------
381	if(co2cond)then
382
383	do ig=1,ngrid
384	if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015
385	DO nw=1,L_NSPECTV
386	albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
387	ENDDO
388	endif
389	enddo ! ngrid
390
391	endif ! co2cond
392
393
394	do ig=1,ngrid ! We calculate here the tracer tendencies. Don't forget that we have to retrieve the dqsurf tendencies we added at the beginning of the routine !
395	dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep - dqsurf(ig,iliq)
396	dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep - dqsurf(ig,iice)
397	enddo
398
399	if (activerunoff) then
400	call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
401	endif
402
403	return
404	end subroutine hydrol

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

Download in other formats:

Original Format