Context Navigation

hydrol.F90 @ 3638

Last change on this file since 3638 was 3291, checked in by mturbet, 20 months ago
improved parameterization of sea ice albedo vs thickness
Property svn:executable set to ``*
File size: 15.9 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 mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
10	use mod_phys_lmdz_para, only : is_master, gather, scatter
11	use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
12	USE surfdat_h
13	use comdiurn_h
14	USE geometry_mod, only: cell_area
15	USE tracer_h
16	! use slab_ice_h
17	USE ocean_slab_mod, ONLY: alb_ice_min,h_alb_ice,snow_min
18	use callkeys_mod, only: albedosnow,alb_ocean,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
19	use radinc_h, only : L_NSPECTV
20
21	implicit none
22
23	!==================================================================
24	!
25	! Purpose
26	! -------
27	! Calculate the surface hydrology and albedo changes.
28	! Both for oceanic and continental regions
29	!
30	! Authors
31	! -------
32	! Adapted from LMDTERRE by B. Charnay (2010). Further
33	! Modifications by R. Wordsworth (2010).
34	! Spectral albedo by M. Turbet (2015).
35	!
36	! Called by
37	! ---------
38	! physiq.F
39	!
40	! Calls
41	! -----
42	! none
43	!
44	! Notes
45	! -----
46	! rnat is terrain type: 0-ocean; 1-continent
47	!
48	!==================================================================
49
50	integer ngrid,nq
51
52	! Inputs
53	! ------
54	real snowlayer
55	parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
56
57	real oceantime ! this is a relaxation timescale for the oceanbulkavg parameterization
58	parameter (oceantime=10243600)
59	logical,save :: oceanbulkavg ! simple parameterization to relax ocean temperatures to the global mean value (crude, 0th order parameterization to mimick ocean heat transport)
60
61	logical,save :: activerunoff ! enable simple runoff scheme?
62	logical,save :: oceanalbvary ! simple parameterization to account for the effect of solar zenith angle on the ocean albedo (for the moment it is not used, but to be included in the future)
63	!$OMP THREADPRIVATE(oceanbulkavg,activerunoff,oceanalbvary)
64
65	! Arguments
66	! ---------
67	real rnat(ngrid) ! rnat is terrain type: 0-ocean; 1-continent
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 albedo(ngrid,L_NSPECTV)
79	real albedo_bareground(ngrid)
80	real albedo_snow_SPECTV(L_NSPECTV)
81	real albedo_co2_ice_SPECTV(L_NSPECTV)
82	real pctsrf_sic(ngrid), sea_ice(ngrid)
83
84	real oceanarea2
85
86	! Output
87	! ------
88	real dqs_hyd(ngrid,nq)
89	real pdtsurf_hyd(ngrid)
90
91	! Local
92	! -----
93	real a,b,E
94	integer ig,iq, nw
95	real fsnoi, subli, fauxo
96	real twater(ngrid)
97	real pcapcal(ngrid)
98	real hicebis(ngrid)
99	real zqsurf(ngrid,nq)
100	real ztsurf(ngrid)
101	real albedo_sic, alb_ice
102	real frac_snow
103
104	integer, save :: ivap, iliq, iice
105	!$OMP THREADPRIVATE(ivap,iliq,iice)
106
107	logical, save :: firstcall=.true.
108	!$OMP THREADPRIVATE(firstcall)
109
110	real :: runoffamount(ngrid)
111	!#ifdef CPP_PARA
112	real :: runoffamount_glo(klon_glo)
113	real :: zqsurf_iliq_glo(klon_glo)
114	real :: rnat_glo(klon_glo)
115	real :: oceanarea_glo
116	real :: cell_area_glo(klon_glo)
117	!#else
118	! real :: runoffamount_glo(ngrid)
119	! real :: zqsurf_iliq_glo(ngrid)
120	!#endif
121
122
123	if(firstcall)then
124
125	oceanbulkavg=.false.
126	oceanalbvary=.false.
127	write(,)"Activate runnoff into oceans?"
128	activerunoff=.false.
129	call getin_p("activerunoff",activerunoff)
130	write(,)" activerunoff = ",activerunoff
131
132	if (activerunoff) then
133	ALLOCATE(runoff(ngrid))
134	runoff(1:ngrid)=0
135	endif
136
137	ivap=igcm_h2o_vap
138	iliq=igcm_h2o_vap
139	iice=igcm_h2o_ice
140
141	write(,) "hydrol: ivap=",ivap
142	write(,) " iliq=",iliq
143	write(,) " iice=",iice
144
145	! Here's the deal: iice is used in place of igcm_h2o_ice both on the
146	! surface and in the atmosphere. ivap is used in
147	! place of igcm_h2o_vap ONLY in the atmosphere, while
148	! iliq is used in place of igcm_h2o_vap ONLY on the
149	! surface.
150	! Soon to be extended to the entire water cycle...
151
152	! LOCAL ocean surface area
153	oceanarea=0.
154	do ig=1,ngrid
155	if(nint(rnat(ig)).eq.0)then
156	oceanarea=oceanarea+cell_area(ig)
157	endif
158	enddo
159
160	if(oceanbulkavg.and.(oceanarea.le.0.))then
161	print*,'How are we supposed to average the ocean'
162	print*,'temperature, when there are no oceans?'
163	call abort
164	endif
165
166	if(activerunoff.and.(oceanarea.le.0.))then
167	print*,'You have enabled runoff, but you have no oceans.'
168	print*,'Where did you think the water was going to go?'
169	call abort
170	endif
171
172	firstcall = .false.
173	endif
174
175	! write (,) "oceanarea", oceanarea
176
177	! add physical tendencies already calculated
178	! ------------------------------------------
179
180	do ig=1,ngrid
181	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
182	pdtsurf_hyd(ig)=0.0
183	do iq=1,nq
184	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
185	enddo
186	enddo
187
188	do ig=1,ngrid
189	do iq=1,nq
190	dqs_hyd(ig,iq) = 0.0
191	enddo
192	enddo
193
194	do ig = 1, ngrid
195
196	! Ocean regions (rnat = 0)
197	! -----------------------
198	if(nint(rnat(ig)).eq.0)then
199
200	! Parameterization (not used for the moment) to compute the effect of solar zenith angle on the albedo
201	! --------------------------
202	!
203	! if(diurnal.and.oceanalbvary)then
204	! fauxo = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
205	! albedo(ig) = 1.1( .03 + .630/( 1. + fauxofauxo))
206	! albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
207	! else
208	!
209	! do nw=1,L_NSPECTV
210	! albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in inifis_mod.F90. Later we could introduce spectral dependency for alb_ocean.
211	! enddo
212	! end if
213
214	! we first start by fixing the albedo of oceanic grid to that of the ocean
215	do nw=1,L_NSPECTV
216	albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in inifis_mod.F90. Later we could introduce spectral dependency for alb_ocean.
217	enddo
218
219	if(ok_slab_ocean) then ! if ocean heat transport param activated
220
221	frac_snow = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/snow_min)) ! Critical snow height (in kg/m2) from ocean_slab_ice routine.
222	! Standard value should be 15kg/m2 (i.e. about 5 cm). Note that in the previous ocean param. (from BC2014), this value was 45kg/m2 (i.e. about 15cm).
223
224	! Albedo final calculation :
225	do nw=1,L_NSPECTV
226	alb_ice=albedo_snow_SPECTV(nw)-(albedo_snow_SPECTV(nw)-alb_ice_min)*exp(-sea_ice(ig)/h_alb_ice) ! this replaces the formulation from BC2014
227	! More details on the parameterization of sea ice albedo vs thickness is provided in the wiki :
228	! https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Slab_ocean_model
229	! sea_ice is the ice thickness (calculated in ocean_slab routine) in kg/m2 ; h_alb_ice is fixed to 275.1kg/m2 i.e. 30cm based on comparisons with Brandt et al. 2005
230	albedo(ig,nw) = pctsrf_sic(ig)* &
231	(albedo_snow_SPECTV(nw)frac_snow + alb_ice(1.0-frac_snow)) &
232	+ (1.-pctsrf_sic(ig))*alb_ocean
233	enddo
234
235	! Oceanic ice height, just for diagnostics
236	hice(ig) = MIN(10.,sea_ice(ig)/rhowater)
237	else !ok_slab_ocean ; here this is the case where we are dealing with a static ocean
238
239
240	! calculate oceanic ice height including the latent heat of ice formation
241	! hice is the height of oceanic ice with a maximum of maxicethick.
242	hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
243	twater(ig) = ztsurf(ig) - hice(ig)RLFTTrhowater/pcapcal(ig) ! this is the temperature water would have if we melted the entire ocean ice layer
244	hicebis(ig) = hice(ig)
245	hice(ig) = 0.
246
247	if(twater(ig) .lt. T_h2O_ice_liq)then
248	E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))pcapcal(ig),RLFTTrhowater*maxicethick)
249	hice(ig) = E/(RLFTT*rhowater)
250	hice(ig) = max(hice(ig),0.0)
251	hice(ig) = min(hice(ig),maxicethick)
252	pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))RLFTTrhowater/pcapcal(ig)/ptimestep
253	do nw=1,L_NSPECTV
254	albedo(ig,nw) = albedo_snow_SPECTV(nw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
255	enddo
256
257	! if (zqsurf(ig,iice).ge.snowlayer) then
258	! albedo(ig) = albedoice
259	! else
260	! albedo(ig) = albedoocean &
261	! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
262	! endif
263
264	else
265
266	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
267	DO nw=1,L_NSPECTV
268	albedo(ig,nw) = alb_ocean
269	ENDDO
270
271	endif
272
273	zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
274	zqsurf(ig,iice) = hice(ig)*rhowater
275
276	endif!(ok_slab_ocean)
277
278
279	! Continental regions (rnat = 1)
280	! -----------------------
281	elseif (nint(rnat(ig)).eq.1) then
282
283	! melt the snow
284	if(ztsurf(ig).gt.T_h2O_ice_liq)then
285	if(zqsurf(ig,iice).gt.1.0e-8)then
286
287	a = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
288	b = zqsurf(ig,iice)
289	fsnoi = min(a,b)
290
291	zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
292	zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi
293
294	! thermal effects
295	pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep
296
297	endif
298	else
299
300	! freeze the water
301	if(zqsurf(ig,iliq).gt.1.0e-8)then
302
303	a = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
304	b = zqsurf(ig,iliq)
305
306	fsnoi = min(a,b)
307
308	zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
309	zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi
310
311	! thermal effects
312	pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep
313
314	endif
315	endif
316
317	! add runoff (to simulate transport of water from continental regions to oceanic regions ; in practice, this prevents liquid water to build up too much on continental regions)
318	if(activerunoff)then
319	runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
320	if(ngrid.gt.1)then ! runoff only exists in 3D
321	if(runoff(ig).ne.0.0)then
322	zqsurf(ig,iliq) = mx_eau_sol
323	! note: runoff is added to ocean at end
324	endif
325	end if
326
327	endif
328
329	! re-calculate continental albedo
330	DO nw=1,L_NSPECTV
331	albedo(ig,nw) = albedo_bareground(ig)
332	ENDDO
333	if (zqsurf(ig,iice).ge.snowlayer) then
334	DO nw=1,L_NSPECTV
335	albedo(ig,nw) = albedo_snow_SPECTV(nw)
336	ENDDO
337	else
338	DO nw=1,L_NSPECTV
339	albedo(ig,nw) = albedo_bareground(ig) &
340	+ (albedo_snow_SPECTV(nw) - albedo_bareground(ig)) &
341	*zqsurf(ig,iice)/snowlayer
342	ENDDO
343	endif
344
345	else
346
347	print*,'Surface type not recognised in hydrol.F!'
348	print*,'Exiting...'
349	call abort
350
351	endif
352
353	end do ! ig=1,ngrid
354
355
356	! simple parameterization to perform crude bulk averaging of temperature in ocean
357	! ----------------------------------------------------
358	if(oceanbulkavg)then
359
360	oceanarea2=0.
361	DO ig=1,ngrid
362	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
363	oceanarea2=oceanarea2+cell_area(ig)*pcapcal(ig)
364	end if
365	END DO
366
367	tsea=0.
368	DO ig=1,ngrid
369	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
370	tsea=tsea+ztsurf(ig)cell_area(ig)pcapcal(ig)/oceanarea2
371	end if
372	END DO
373
374	DO ig=1,ngrid
375	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0))then
376	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
377	end if
378	END DO
379
380	print*,'Mean ocean temperature = ',tsea,' K'
381
382	endif
383
384	! shove all the runoff water into the ocean
385	! -----------------------------------------
386	if(activerunoff)then
387
388	! totalrunoff=0.
389	do ig=1,ngrid
390	runoffamount(ig) = cell_area(ig)*runoff(ig)
391	! if (nint(rnat(ig)).eq.1) then
392	! totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
393	! endif
394	enddo
395
396	! collect on the full grid
397	call gather(runoffamount,runoffamount_glo)
398	call gather(zqsurf(1:ngrid,iliq),zqsurf_iliq_glo)
399	call gather(rnat,rnat_glo)
400	call gather(cell_area,cell_area_glo)
401
402	if (is_master) then
403	totalrunoff=0.
404	oceanarea_glo=0.
405	do ig=1,klon_glo
406	if (nint(rnat_glo(ig)).eq.1) then
407	totalrunoff = totalrunoff + runoffamount_glo(ig)
408	endif
409	if (nint(rnat_glo(ig)).eq.0) then
410	oceanarea_glo = oceanarea_glo + cell_area_glo(ig)
411	endif
412	enddo
413
414	do ig=1,klon_glo
415	if (nint(rnat_glo(ig)).eq.0) then
416	zqsurf_iliq_glo(ig) = zqsurf_iliq_glo(ig) + &
417	totalrunoff/oceanarea_glo
418	endif
419	enddo
420
421	endif! is_master
422
423	! scatter the field back on all processes
424	call scatter(zqsurf_iliq_glo,zqsurf(1:ngrid,iliq))
425
426
427
428	! do ig=1,ngrid
429	! if (nint(rnat(ig)).eq.0) then
430	! zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
431	! totalrunoff/oceanarea
432	! endif
433	! enddo
434
435	endif !activerunoff
436
437	! Re-add the albedo effects of CO2 ice if necessary
438	! -------------------------------------------------
439	if(co2cond)then
440
441	do ig=1,ngrid
442	if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015
443	DO nw=1,L_NSPECTV
444	albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
445	ENDDO
446	endif
447	enddo ! ngrid
448
449	endif ! co2cond
450
451
452	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 !
453	dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep - dqsurf(ig,iliq)
454	dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep - dqsurf(ig,iice)
455	enddo
456
457	if (activerunoff) then
458	call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
459	endif
460
461	return
462	end subroutine hydrol

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

Download in other formats:

Original Format