Context Navigation

hydrol.F90 @ 1243

Last change on this file since 1243 was 965, checked in by emillour, 12 years ago

Common dynamics and generic/universal GCM:

LMDZ.COMMON: minor bug fix on the computation of physics mesh area in gcm.F
LMDZ.UNIVERSAL: missing clean initialization of tab_cntrl(:) array in phyredem.F90
LMDZ.GENERIC: minor bug fix in hydrol.F90, only output runoff if it is used. Update output routines so that all outputs files (stats, diagfi.nc, diagsoil.nc, diagspecIR.nc and diagspecVI.nc) can be generated when running LMDZ.UNIVERSAL in MPI mode.

Property svn:executable set to *

File size: 10.4 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
[965]	112	if (activerunoff) ALLOCATE(runoff(ngrid))
[787]	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
[965]	362	if (activerunoff) then
	363	call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
	364	endif
[253]	365
	366	return
	367	end subroutine hydrol

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

Download in other formats:

Original Format