Context Navigation

hydrol.F90 @ 2921

Last change on this file since 2921 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

Rev	Line
[1482]	1	subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf, &
[253]	2	qsurf,dqsurf,dqs_hyd,pcapcal, &
[1482]	3	albedo,albedo_bareground, &
	4	albedo_snow_SPECTV,albedo_co2_ice_SPECTV, &
	5	mu0,pdtsurf,pdtsurf_hyd,hice, &
[1297]	6	pctsrf_sic,sea_ice)
[253]	7
[1521]	8	use ioipsl_getin_p_mod, only: getin_p
[650]	9	use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
[787]	10	USE surfdat_h
	11	use comdiurn_h
[1543]	12	USE geometry_mod, only: cell_area
[787]	13	USE tracer_h
[1297]	14	use slab_ice_h
[2482]	15	use callkeys_mod, only: albedosnow,alb_ocean,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
[1482]	16	use radinc_h, only : L_NSPECTV
[253]	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
[1482]	29	! Modifications by R. Wordsworth (2010).
	30	! Spectral albedo by M. Turbet (2015).
[253]	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
[787]	46	integer ngrid,nq
	47
[253]	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
[305]	53	parameter (oceantime=10243600)
[253]	54
[875]	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?
[1315]	58	!$OMP THREADPRIVATE(oceanbulkavg,activerunoff,oceanalbvary)
[253]	59
	60	! Arguments
	61	! ---------
[1297]	62	real rnat(ngrid) ! I changed this to integer (RW)
[787]	63	real,dimension(:),allocatable,save :: runoff
[253]	64	real totalrunoff, tsea, oceanarea
	65	save oceanarea
[1315]	66	!$OMP THREADPRIVATE(runoff,oceanarea)
[253]	67
	68	real ptimestep
[787]	69	real mu0(ngrid)
	70	real qsurf(ngrid,nq), tsurf(ngrid)
	71	real dqsurf(ngrid,nq), pdtsurf(ngrid)
	72	real hice(ngrid)
[1482]	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)
[1297]	77	real pctsrf_sic(ngrid), sea_ice(ngrid)
[253]	78
	79	real oceanarea2
	80
	81	! Output
	82	! ------
[787]	83	real dqs_hyd(ngrid,nq)
	84	real pdtsurf_hyd(ngrid)
[253]	85
	86	! Local
	87	! -----
	88	real a,b,E
[1482]	89	integer ig,iq, nw
[253]	90	real fsnoi, subli, fauxo
[787]	91	real twater(ngrid)
	92	real pcapcal(ngrid)
	93	real hicebis(ngrid)
	94	real zqsurf(ngrid,nq)
	95	real ztsurf(ngrid)
[1297]	96	real albedo_sic, alb_ice
	97	real zfra
[253]	98
[863]	99	integer, save :: ivap, iliq, iice
[1315]	100	!$OMP THREADPRIVATE(ivap,iliq,iice)
[253]	101
[863]	102	logical, save :: firstcall
[1315]	103	!$OMP THREADPRIVATE(firstcall)
[253]	104
	105	data firstcall /.true./
	106
	107
	108	if(firstcall)then
	109
[875]	110	oceanbulkavg=.false.
	111	oceanalbvary=.false.
	112	write(,)"Activate runnoff into oceans?"
	113	activerunoff=.false.
[1315]	114	call getin_p("activerunoff",activerunoff)
[875]	115	write(,)" activerunoff = ",activerunoff
	116
	117
	118
[1537]	119	if (activerunoff) then
	120	ALLOCATE(runoff(ngrid))
	121	runoff(1:ngrid)=0
	122	endif
[787]	123
[253]	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.
[787]	141	do ig=1,ngrid
[1297]	142	if(nint(rnat(ig)).eq.0)then
[1542]	143	oceanarea=oceanarea+cell_area(ig)
[253]	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
[787]	165	do ig=1,ngrid
[253]	166	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
	167	pdtsurf_hyd(ig)=0.0
[787]	168	do iq=1,nq
[253]	169	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
	170	enddo
	171	enddo
	172
[787]	173	do ig=1,ngrid
	174	do iq=1,nq
[253]	175	dqs_hyd(ig,iq) = 0.0
	176	enddo
	177	enddo
	178
[787]	179	do ig = 1, ngrid
[253]	180
	181	! Ocean
	182	! -----
[1297]	183	if(nint(rnat(ig)).eq.0)then
[253]	184
	185	! re-calculate oceanic albedo
[1297]	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
[1482]	191	do nw=1,L_NSPECTV
[2482]	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.
[1482]	193	enddo
[1297]	194	! end if
[253]	195
[1297]	196
[1482]	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
[1297]	208
[1482]	209	! Oceanic ice height, just for diagnostics
	210	hice(ig) = MIN(10.,sea_ice(ig)/rhowater)
	211	else !ok_slab_ocean
[1297]	212
	213
[253]	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.
[1482]	216	hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
[253]	217
[1482]	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.
[253]	223
[1482]	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
[253]	233
[1482]	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
[253]	240
[1482]	241	else
[253]	242
[1482]	243	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
	244	DO nw=1,L_NSPECTV
	245	albedo(ig,nw) = alb_ocean
	246	ENDDO
[253]	247
[1482]	248	endif
[253]	249
[1482]	250	zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
	251	zqsurf(ig,iice) = hice(ig)*rhowater
[253]	252
[1482]	253	endif!(ok_slab_ocean)
[253]	254
[1297]	255
[253]	256	! Continent
	257	! ---------
[1297]	258	elseif (nint(rnat(ig)).eq.1) then
[253]	259
	260	! melt the snow
[650]	261	if(ztsurf(ig).gt.T_h2O_ice_liq)then
[253]	262	if(zqsurf(ig,iice).gt.1.0e-8)then
	263
[650]	264	a = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
[253]	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
[650]	280	a = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
[253]	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)
[787]	298	if(ngrid.gt.1)then ! runoff only exists in 3D
[253]	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
[1482]	308	DO nw=1,L_NSPECTV
	309	albedo(ig,nw) = albedo_bareground(ig)
	310	ENDDO
[253]	311	if (zqsurf(ig,iice).ge.snowlayer) then
[1482]	312	DO nw=1,L_NSPECTV
	313	albedo(ig,nw) = albedo_snow_SPECTV(nw)
	314	ENDDO
[253]	315	else
[1482]	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
[253]	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
[787]	331	end do ! ig=1,ngrid
[253]	332
	333	! perform crude bulk averaging of temperature in ocean
	334	! ----------------------------------------------------
	335	if(oceanbulkavg)then
	336
	337	oceanarea2=0.
[787]	338	DO ig=1,ngrid
[1297]	339	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
[1542]	340	oceanarea2=oceanarea2+cell_area(ig)*pcapcal(ig)
[253]	341	end if
	342	END DO
	343
	344	tsea=0.
[787]	345	DO ig=1,ngrid
[1297]	346	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
[1542]	347	tsea=tsea+ztsurf(ig)cell_area(ig)pcapcal(ig)/oceanarea2
[253]	348	end if
	349	END DO
	350
[787]	351	DO ig=1,ngrid
[1297]	352	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0))then
[253]	353	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
	354	end if
	355	END DO
	356
[305]	357	print*,'Mean ocean temperature = ',tsea,' K'
[253]	358
	359	endif
	360
	361	! shove all the runoff water into the ocean
	362	! -----------------------------------------
	363	if(activerunoff)then
	364
	365	totalrunoff=0.
[787]	366	do ig=1,ngrid
[1297]	367	if (nint(rnat(ig)).eq.1) then
[1542]	368	totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
[253]	369	endif
	370	enddo
	371
[787]	372	do ig=1,ngrid
[1297]	373	if (nint(rnat(ig)).eq.0) then
[253]	374	zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
	375	totalrunoff/oceanarea
	376	endif
	377	enddo
	378
[863]	379	endif
[253]	380
[863]	381
[253]	382	! Re-add the albedo effects of CO2 ice if necessary
	383	! -------------------------------------------------
	384	if(co2cond)then
	385
[787]	386	do ig=1,ngrid
[1482]	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
[253]	391	endif
[1482]	392	enddo ! ngrid
[253]	393
[1482]	394	endif ! co2cond
[253]	395
	396
[1484]	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)
[253]	400	enddo
	401
[965]	402	if (activerunoff) then
	403	call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
	404	endif
[253]	405
	406	return
	407	end subroutine hydrol

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

Download in other formats:

Original Format