Context Navigation

hydrol.F90 @ 1330

Last change on this file since 1330 was 1315, checked in by milmd, 11 years ago
LMDZ.GENERIC. OpenMP directives added in generic physic. When running in pure OpenMP or hybrid OpenMP/MPI, may have some bugs with condense_cloud and wstats routines.
Property svn:executable set to ``*
File size: 11.3 KB

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

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

Download in other formats:

Original Format