Context Navigation

hydrol.F90 @ 1579

Last change on this file since 1579 was 1543, checked in by emillour, 10 years ago

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

dyn3d:
adapted gcm.F so that all physics initializations are now done in iniphysiq.

dyn3dpar:
adapted gcm.F so that all physics initializations are now done in iniphysiq.
updated calfis_p.F to follow up with changes.
copied over updated "bands.F90" from LMDZ5.

dynphy_lonlat:
calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines

phy_common:
added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
added "physics_distribution_mod.F90"
updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.

phy[venus/titan/mars/std]:
removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

Property svn:executable set to *

File size: 12.6 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 geometry_mod, only: cell_area
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) then
120	ALLOCATE(runoff(ngrid))
121	runoff(1:ngrid)=0
122	endif
123
124	ivap=igcm_h2o_vap
125	iliq=igcm_h2o_vap
126	iice=igcm_h2o_ice
127
128	write(,) "hydrol: ivap=",ivap
129	write(,) " iliq=",iliq
130	write(,) " iice=",iice
131
132	! Here's the deal: iice is used in place of igcm_h2o_ice both on the
133	! surface and in the atmosphere. ivap is used in
134	! place of igcm_h2o_vap ONLY in the atmosphere, while
135	! iliq is used in place of igcm_h2o_vap ONLY on the
136	! surface.
137	! Soon to be extended to the entire water cycle...
138
139	! Total ocean surface area
140	oceanarea=0.
141	do ig=1,ngrid
142	if(nint(rnat(ig)).eq.0)then
143	oceanarea=oceanarea+cell_area(ig)
144	endif
145	enddo
146
147	if(oceanbulkavg.and.(oceanarea.le.0.))then
148	print*,'How are we supposed to average the ocean'
149	print*,'temperature, when there are no oceans?'
150	call abort
151	endif
152
153	if(activerunoff.and.(oceanarea.le.0.))then
154	print*,'You have enabled runoff, but you have no oceans.'
155	print*,'Where did you think the water was going to go?'
156	call abort
157	endif
158
159	firstcall = .false.
160	endif
161
162	! add physical tendencies already calculated
163	! ------------------------------------------
164
165	do ig=1,ngrid
166	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
167	pdtsurf_hyd(ig)=0.0
168	do iq=1,nq
169	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
170	enddo
171	enddo
172
173	do ig=1,ngrid
174	do iq=1,nq
175	dqs_hyd(ig,iq) = 0.0
176	enddo
177	enddo
178
179	do ig = 1, ngrid
180
181	! Ocean
182	! -----
183	if(nint(rnat(ig)).eq.0)then
184
185	! re-calculate oceanic albedo
186	! if(diurnal.and.oceanalbvary)then
187	! fauxo = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
188	! albedo(ig) = 1.1( .03 + .630/( 1. + fauxofauxo))
189	! albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
190	! else
191	do nw=1,L_NSPECTV
192	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.
193	enddo
194	! end if
195
196
197	if(ok_slab_ocean) then
198
199	zfra = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/45.0)) ! Snow Fraction (Critical height 45kg/m2~15cm)
200	alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) & ! Ice Albedo
201	*exp(-sea_ice(ig)/h_alb_ice)
202	! Albedo final calculation :
203	do nw=1,L_NSPECTV
204	albedo(ig,nw) = pctsrf_sic(ig)* &
205	(albedo_snow_SPECTV(nw)zfra + alb_ice(1.0-zfra)) &
206	+ (1.-pctsrf_sic(ig))*alb_ocean
207	enddo
208
209	! Oceanic ice height, just for diagnostics
210	hice(ig) = MIN(10.,sea_ice(ig)/rhowater)
211	else !ok_slab_ocean
212
213
214	! calculate oceanic ice height including the latent heat of ice formation
215	! hice is the height of oceanic ice with a maximum of maxicethick.
216	hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
217
218	! twater(ig) = tsurf(ig) + ptimestep*zdtsurf(ig) &
219	twater(ig) = ztsurf(ig) - hice(ig)RLFTTrhowater/pcapcal(ig)
220	! this is temperature water would have if we melted entire ocean ice layer
221	hicebis(ig) = hice(ig)
222	hice(ig) = 0.
223
224	if(twater(ig) .lt. T_h2O_ice_liq)then
225	E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))pcapcal(ig),RLFTTrhowater*maxicethick)
226	hice(ig) = E/(RLFTT*rhowater)
227	hice(ig) = max(hice(ig),0.0)
228	hice(ig) = min(hice(ig),maxicethick)
229	pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))RLFTTrhowater/pcapcal(ig)/ptimestep
230	do nw=1,L_NSPECTV
231	albedo(ig,nw) = albedo_snow_SPECTV(nw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
232	enddo
233
234	! if (zqsurf(ig,iice).ge.snowlayer) then
235	! albedo(ig) = albedoice
236	! else
237	! albedo(ig) = albedoocean &
238	! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
239	! endif
240
241	else
242
243	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
244	DO nw=1,L_NSPECTV
245	albedo(ig,nw) = alb_ocean
246	ENDDO
247
248	endif
249
250	zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
251	zqsurf(ig,iice) = hice(ig)*rhowater
252
253	endif!(ok_slab_ocean)
254
255
256	! Continent
257	! ---------
258	elseif (nint(rnat(ig)).eq.1) then
259
260	! melt the snow
261	if(ztsurf(ig).gt.T_h2O_ice_liq)then
262	if(zqsurf(ig,iice).gt.1.0e-8)then
263
264	a = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
265	b = zqsurf(ig,iice)
266	fsnoi = min(a,b)
267
268	zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
269	zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi
270
271	! thermal effects
272	pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep
273
274	endif
275	else
276
277	! freeze the water
278	if(zqsurf(ig,iliq).gt.1.0e-8)then
279
280	a = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
281	b = zqsurf(ig,iliq)
282
283	fsnoi = min(a,b)
284
285	zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
286	zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi
287
288	! thermal effects
289	pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep
290
291	endif
292	endif
293
294	! deal with runoff
295	if(activerunoff)then
296
297	runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
298	if(ngrid.gt.1)then ! runoff only exists in 3D
299	if(runoff(ig).ne.0.0)then
300	zqsurf(ig,iliq) = mx_eau_sol
301	! runoff is added to ocean at end
302	endif
303	end if
304
305	endif
306
307	! re-calculate continental albedo
308	DO nw=1,L_NSPECTV
309	albedo(ig,nw) = albedo_bareground(ig)
310	ENDDO
311	if (zqsurf(ig,iice).ge.snowlayer) then
312	DO nw=1,L_NSPECTV
313	albedo(ig,nw) = albedo_snow_SPECTV(nw)
314	ENDDO
315	else
316	DO nw=1,L_NSPECTV
317	albedo(ig,nw) = albedo_bareground(ig) &
318	+ (albedo_snow_SPECTV(nw) - albedo_bareground(ig)) &
319	*zqsurf(ig,iice)/snowlayer
320	ENDDO
321	endif
322
323	else
324
325	print*,'Surface type not recognised in hydrol.F!'
326	print*,'Exiting...'
327	call abort
328
329	endif
330
331	end do ! ig=1,ngrid
332
333	! perform crude bulk averaging of temperature in ocean
334	! ----------------------------------------------------
335	if(oceanbulkavg)then
336
337	oceanarea2=0.
338	DO ig=1,ngrid
339	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
340	oceanarea2=oceanarea2+cell_area(ig)*pcapcal(ig)
341	end if
342	END DO
343
344	tsea=0.
345	DO ig=1,ngrid
346	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
347	tsea=tsea+ztsurf(ig)cell_area(ig)pcapcal(ig)/oceanarea2
348	end if
349	END DO
350
351	DO ig=1,ngrid
352	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0))then
353	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
354	end if
355	END DO
356
357	print*,'Mean ocean temperature = ',tsea,' K'
358
359	endif
360
361	! shove all the runoff water into the ocean
362	! -----------------------------------------
363	if(activerunoff)then
364
365	totalrunoff=0.
366	do ig=1,ngrid
367	if (nint(rnat(ig)).eq.1) then
368	totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
369	endif
370	enddo
371
372	do ig=1,ngrid
373	if (nint(rnat(ig)).eq.0) then
374	zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
375	totalrunoff/oceanarea
376	endif
377	enddo
378
379	endif
380
381
382	! Re-add the albedo effects of CO2 ice if necessary
383	! -------------------------------------------------
384	if(co2cond)then
385
386	do ig=1,ngrid
387	if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015
388	DO nw=1,L_NSPECTV
389	albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
390	ENDDO
391	endif
392	enddo ! ngrid
393
394	endif ! co2cond
395
396
397	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 !
398	dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep - dqsurf(ig,iliq)
399	dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep - dqsurf(ig,iice)
400	enddo
401
402	if (activerunoff) then
403	call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
404	endif
405
406	return
407	end subroutine hydrol

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

Download in other formats:

Original Format