Context Navigation

hydrol.F90 @ 2613

Last change on this file since 2613 was 2482, checked in by yjaziri, 4 years ago

Generic GCM:
Clean convadj.F90 specific CO2 Mars convection
Add alb_ocean used in hydrol.F90 as option in .def files
Add kmixmin 1D minimum eddy mix coeff for turbdiff as rcm1d.def option
and comment lines to help coding specific eddy mix coeff in turbdiff with Earth example

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,alb_ocean,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 inifis_mod.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