Context Navigation

hydrol.F90 @ 3267

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

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

Download in other formats:

Original Format