Context Navigation

hydrol.F90 @ 1351

Last change on this file since 1351 was 1315, checked in by milmd, 10 years ago
LMDZ.GENERIC. OpenMP directives added in generic physic. When running in pure OpenMP or hybrid OpenMP/MPI, may have some bugs with condense_cloud and wstats routines.
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:

Original Format