Context Navigation

hydrol.F90 @ 918

Last change on this file since 918 was 875, checked in by jleconte, 12 years ago

11/02/2013 == JL

Updated moist convection scheme to handle situations with a large water vapor content
Added a keyword to enable ocean runoff in callphys.def (activerunoff)

Property svn:executable set to *

File size: 10.3 KB

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

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

Download in other formats:

Original Format