Context Navigation

hydrol.F90 @ 537

Last change on this file since 537 was 368, checked in by rwordsworth, 13 years ago
.nc saves in physiq.F90 adapted to new diagfi.def format. More info added to gases.def warnings.
Property svn:executable set to ``*
File size: 10.1 KB

Rev	Line
[253]	1	subroutine hydrol(ptimestep,rnat,tsurf, &
	2	qsurf,dqsurf,dqs_hyd,pcapcal, &
	3	albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
	4
	5	use watercommon_h, only: To, RLFTT, rhowater, mx_eau_sol
	6
	7	implicit none
	8
	9	!==================================================================
	10	!
	11	! Purpose
	12	! -------
	13	! Calculate the surface hydrology and albedo changes.
	14	!
	15	! Authors
	16	! -------
	17	! Adapted from LMDTERRE by B. Charnay (2010). Further
	18	! modifications by R. Wordsworth (2010).
	19	!
	20	! Called by
	21	! ---------
	22	! physiq.F
	23	!
	24	! Calls
	25	! -----
	26	! none
	27	!
	28	! Notes
	29	! -----
	30	! rnat is terrain type: 0-ocean; 1-continent
	31	!
	32	!==================================================================
	33
	34	#include "dimensions.h"
	35	#include "dimphys.h"
	36	#include "comcstfi.h"
	37	#include "callkeys.h"
	38	#include "tracer.h"
	39	#include "fisice.h"
	40	#include "comgeomfi.h"
	41	#include "comdiurn.h"
	42	#include "surfdat.h"
	43
	44	! Inputs
	45	! ------
	46	real albedoocean
	47	parameter (albedoocean=0.07)
	48	real albedoice
	49	save albedoice
	50
	51	real snowlayer
	52	parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
	53	real oceantime
[305]	54	parameter (oceantime=10243600)
[253]	55
	56	logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
	57	logical activerunoff ! enable simple runoff scheme?
	58	logical oceanalbvary ! ocean albedo varies with the diurnal cycle?
	59
	60	! Arguments
	61	! ---------
	62	integer rnat(ngridmx) ! I changed this to integer (RW)
	63	real runoff(ngridmx)
	64	save runoff
	65	real totalrunoff, tsea, oceanarea
	66	save oceanarea
	67
	68	real ptimestep
	69	real mu0(ngridmx)
	70	real qsurf(ngridmx,nqmx), tsurf(ngridmx)
	71	real dqsurf(ngridmx,nqmx), pdtsurf(ngridmx)
	72	real hice(ngridmx)
	73	real albedo0(ngridmx), albedo(ngridmx)
	74
	75	real oceanarea2
	76
	77	! Output
	78	! ------
	79	real dqs_hyd(ngridmx,nqmx)
	80	real pdtsurf_hyd(ngridmx)
	81
	82	! Local
	83	! -----
	84	real a,b,E
	85	integer ig,iq, icap ! wld like to remove icap
	86	real fsnoi, subli, fauxo
	87	real twater(ngridmx)
	88	real pcapcal(ngridmx)
	89	real hicebis(ngridmx)
	90	real zqsurf(ngridmx,nqmx)
	91	real ztsurf(ngridmx)
	92
	93	integer ivap, iliq, iice
	94	save ivap, iliq, iice
	95
	96	logical firstcall
	97	save firstcall
	98
	99	data firstcall /.true./
	100
	101	oceanbulkavg=.false.
[368]	102	activerunoff=.false.
[253]	103	oceanalbvary=.false.
	104
	105	if(firstcall)then
	106
	107	ivap=igcm_h2o_vap
	108	iliq=igcm_h2o_vap
	109	iice=igcm_h2o_ice
	110
	111	write(,) "hydrol: ivap=",ivap
	112	write(,) " iliq=",iliq
	113	write(,) " iice=",iice
	114
	115	! Here's the deal: iice is used in place of igcm_h2o_ice both on the
	116	! surface and in the atmosphere. ivap is used in
	117	! place of igcm_h2o_vap ONLY in the atmosphere, while
	118	! iliq is used in place of igcm_h2o_vap ONLY on the
	119	! surface.
	120
	121	! Soon to be extended to the entire water cycle...
	122
	123	! Ice albedo = snow albedo for now
	124	albedoice=albedosnow
	125
	126	! Total ocean surface area
	127	oceanarea=0.
	128	do ig=1,ngridmx
	129	if(rnat(ig).eq.0)then
	130	oceanarea=oceanarea+area(ig)
	131	endif
	132	enddo
	133
	134	if(oceanbulkavg.and.(oceanarea.le.0.))then
	135	print*,'How are we supposed to average the ocean'
	136	print*,'temperature, when there are no oceans?'
	137	call abort
	138	endif
	139
	140	if(activerunoff.and.(oceanarea.le.0.))then
	141	print*,'You have enabled runoff, but you have no oceans.'
	142	print*,'Where did you think the water was going to go?'
	143	call abort
	144	endif
	145
	146	firstcall = .false.
	147	endif
	148
	149	! add physical tendencies already calculated
	150	! ------------------------------------------
	151
	152	do ig=1,ngridmx
	153	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
	154	pdtsurf_hyd(ig)=0.0
	155	do iq=1,nqmx
	156	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
	157	enddo
	158	enddo
	159
	160	do ig=1,ngridmx
	161	do iq=1,nqmx
	162	dqs_hyd(ig,iq) = 0.0
	163	enddo
	164	enddo
	165
	166	do ig = 1, ngridmx
	167
	168	! Ocean
	169	! -----
	170	if(rnat(ig).eq.0)then
	171
	172	! re-calculate oceanic albedo
	173	if(diurnal.and.oceanalbvary)then
	174	fauxo = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
	175	albedo(ig) = 1.1( .03 + .630/( 1. + fauxofauxo))
	176	albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
	177	else
	178	albedo(ig) = albedoocean ! modif Benjamin
	179	end if
	180
	181	! calculate oceanic ice height including the latent heat of ice formation
	182	! hice is the height of oceanic ice with a maximum of maxicethick.
	183	hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
	184
	185	! twater(ig) = tsurf(ig) + ptimestep*zdtsurf(ig) &
	186	twater(ig) = ztsurf(ig) - hice(ig)RLFTTrhowater/pcapcal(ig)
	187	! this is temperature water would have if we melted entire ocean ice layer
	188	hicebis(ig) = hice(ig)
	189	hice(ig) = 0.
	190
	191	if(twater(ig) .lt. To)then
	192	E=min((To+Tsaldiff-twater(ig))pcapcal(ig),RLFTTrhowater*maxicethick)
	193	hice(ig) = E/(RLFTT*rhowater)
	194	hice(ig) = max(hice(ig),0.0)
	195	hice(ig) = min(hice(ig),maxicethick)
	196	pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))RLFTT &
	197	rhowater/pcapcal(ig)/ptimestep
	198	albedo(ig) = albedoice
	199
	200	! if (zqsurf(ig,iice).ge.snowlayer) then
	201	! albedo(ig) = albedoice
	202	! else
	203	! albedo(ig) = albedoocean &
	204	! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
	205	! endif
	206
	207	else
	208
	209	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
	210	albedo(ig) = albedoocean
	211
	212	endif
	213
	214	zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
	215	zqsurf(ig,iice) = hice(ig)*rhowater
	216
	217
	218	! Continent
	219	! ---------
	220	elseif (rnat(ig).eq.1) then
	221
	222	! melt the snow
	223	if(ztsurf(ig).gt.To)then
	224	if(zqsurf(ig,iice).gt.1.0e-8)then
	225
	226	a = (ztsurf(ig)-To)*pcapcal(ig)/RLFTT
	227	b = zqsurf(ig,iice)
	228	fsnoi = min(a,b)
	229
	230	zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
	231	zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi
	232
	233	! thermal effects
	234	pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep
	235
	236	endif
	237	else
	238
	239	! freeze the water
	240	if(zqsurf(ig,iliq).gt.1.0e-8)then
	241
	242	a = -(ztsurf(ig)-To)*pcapcal(ig)/RLFTT
	243	b = zqsurf(ig,iliq)
	244
	245	fsnoi = min(a,b)
	246
	247	zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
	248	zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi
	249
	250	! thermal effects
	251	pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep
	252
	253	endif
	254	endif
	255
	256	! deal with runoff
	257	if(activerunoff)then
	258
	259	runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
	260	if(ngridmx.gt.1)then ! runoff only exists in 3D
	261	if(runoff(ig).ne.0.0)then
	262	zqsurf(ig,iliq) = mx_eau_sol
	263	! runoff is added to ocean at end
	264	endif
	265	end if
	266
	267	endif
	268
	269	! re-calculate continental albedo
	270	albedo(ig) = albedo0(ig) ! albedo0 = base values
	271	if (zqsurf(ig,iice).ge.snowlayer) then
	272	albedo(ig) = albedosnow
	273	else
	274	albedo(ig) = albedo0(ig) &
	275	+ (albedosnow - albedo0(ig))*zqsurf(ig,iice)/snowlayer
	276	endif
	277
	278	else
	279
	280	print*,'Surface type not recognised in hydrol.F!'
	281	print*,'Exiting...'
	282	call abort
	283
	284	endif
	285
	286	end do ! ig=1,ngridmx
	287
	288	! perform crude bulk averaging of temperature in ocean
	289	! ----------------------------------------------------
	290	if(oceanbulkavg)then
	291
	292	oceanarea2=0.
	293	DO ig=1,ngridmx
	294	if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
	295	oceanarea2=oceanarea2+area(ig)*pcapcal(ig)
	296	end if
	297	END DO
	298
	299	tsea=0.
	300	DO ig=1,ngridmx
	301	if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
	302	tsea=tsea+ztsurf(ig)area(ig)pcapcal(ig)/oceanarea2
	303	end if
	304	END DO
	305
	306	DO ig=1,ngridmx
	307	if((rnat(ig).eq.0).and.(hice(ig).eq.0))then
	308	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
	309	end if
	310	END DO
	311
[305]	312	print*,'Mean ocean temperature = ',tsea,' K'
[253]	313
	314	endif
	315
	316	! shove all the runoff water into the ocean
	317	! -----------------------------------------
	318	if(activerunoff)then
	319
	320	totalrunoff=0.
	321	do ig=1,ngridmx
	322	if (rnat(ig).eq.1) then
	323	totalrunoff = totalrunoff + area(ig)*runoff(ig)
	324	endif
	325	enddo
	326
	327	do ig=1,ngridmx
	328	if (rnat(ig).eq.0) then
	329	zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
	330	totalrunoff/oceanarea
	331	endif
	332	enddo
	333
	334	endif
	335
	336	! Re-add the albedo effects of CO2 ice if necessary
	337	! -------------------------------------------------
	338	if(co2cond)then
	339
	340	icap=1
	341	do ig=1,ngridmx
	342	if (qsurf(ig,igcm_co2_ice).gt.0) then
	343	albedo(ig) = albedice(icap)
	344	endif
	345	enddo
	346
	347	endif
	348
	349
	350	do ig=1,ngridmx
	351	dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep
	352	dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep
	353	enddo
	354
	355	call writediagfi(ngridmx,'runoff','Runoff amount',' ',2,runoff)
	356
	357	return
	358	end subroutine hydrol

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

Download in other formats:

Original Format