Context Navigation

hydrol.F90 @ 3469

Last change on this file since 3469 was 3291, checked in by mturbet, 9 months ago
improved parameterization of sea ice albedo vs thickness
Property svn:executable set to ``*
File size: 15.9 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
[3100]	9	use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
	10	use mod_phys_lmdz_para, only : is_master, gather, scatter
[650]	11	use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
[787]	12	USE surfdat_h
	13	use comdiurn_h
[1543]	14	USE geometry_mod, only: cell_area
[787]	15	USE tracer_h
[3100]	16	! use slab_ice_h
[3291]	17	USE ocean_slab_mod, ONLY: alb_ice_min,h_alb_ice,snow_min
[2482]	18	use callkeys_mod, only: albedosnow,alb_ocean,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
[1482]	19	use radinc_h, only : L_NSPECTV
[253]	20
	21	implicit none
	22
	23	!==================================================================
	24	!
	25	! Purpose
	26	! -------
	27	! Calculate the surface hydrology and albedo changes.
[3267]	28	! Both for oceanic and continental regions
[253]	29	!
	30	! Authors
	31	! -------
	32	! Adapted from LMDTERRE by B. Charnay (2010). Further
[1482]	33	! Modifications by R. Wordsworth (2010).
	34	! Spectral albedo by M. Turbet (2015).
[253]	35	!
	36	! Called by
	37	! ---------
	38	! physiq.F
	39	!
	40	! Calls
	41	! -----
	42	! none
	43	!
	44	! Notes
	45	! -----
	46	! rnat is terrain type: 0-ocean; 1-continent
	47	!
	48	!==================================================================
	49
[787]	50	integer ngrid,nq
	51
[253]	52	! Inputs
	53	! ------
	54	real snowlayer
	55	parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
[3266]	56
	57	real oceantime ! this is a relaxation timescale for the oceanbulkavg parameterization
[305]	58	parameter (oceantime=10243600)
[3266]	59	logical,save :: oceanbulkavg ! simple parameterization to relax ocean temperatures to the global mean value (crude, 0th order parameterization to mimick ocean heat transport)
	60
[875]	61	logical,save :: activerunoff ! enable simple runoff scheme?
[3266]	62	logical,save :: oceanalbvary ! simple parameterization to account for the effect of solar zenith angle on the ocean albedo (for the moment it is not used, but to be included in the future)
[1315]	63	!$OMP THREADPRIVATE(oceanbulkavg,activerunoff,oceanalbvary)
[253]	64
	65	! Arguments
	66	! ---------
[3266]	67	real rnat(ngrid) ! rnat is terrain type: 0-ocean; 1-continent
[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)
[1482]	78	real albedo(ngrid,L_NSPECTV)
	79	real albedo_bareground(ngrid)
	80	real albedo_snow_SPECTV(L_NSPECTV)
	81	real albedo_co2_ice_SPECTV(L_NSPECTV)
[1297]	82	real pctsrf_sic(ngrid), sea_ice(ngrid)
[253]	83
	84	real oceanarea2
	85
	86	! Output
	87	! ------
[787]	88	real dqs_hyd(ngrid,nq)
	89	real pdtsurf_hyd(ngrid)
[253]	90
	91	! Local
	92	! -----
	93	real a,b,E
[1482]	94	integer ig,iq, nw
[253]	95	real fsnoi, subli, fauxo
[787]	96	real twater(ngrid)
	97	real pcapcal(ngrid)
	98	real hicebis(ngrid)
	99	real zqsurf(ngrid,nq)
	100	real ztsurf(ngrid)
[1297]	101	real albedo_sic, alb_ice
[3268]	102	real frac_snow
[253]	103
[863]	104	integer, save :: ivap, iliq, iice
[1315]	105	!$OMP THREADPRIVATE(ivap,iliq,iice)
[253]	106
[3100]	107	logical, save :: firstcall=.true.
[1315]	108	!$OMP THREADPRIVATE(firstcall)
[253]	109
[3267]	110	real :: runoffamount(ngrid)
[3100]	111	!#ifdef CPP_PARA
	112	real :: runoffamount_glo(klon_glo)
	113	real :: zqsurf_iliq_glo(klon_glo)
	114	real :: rnat_glo(klon_glo)
	115	real :: oceanarea_glo
	116	real :: cell_area_glo(klon_glo)
	117	!#else
	118	! real :: runoffamount_glo(ngrid)
	119	! real :: zqsurf_iliq_glo(ngrid)
	120	!#endif
[253]	121
	122
	123	if(firstcall)then
	124
[875]	125	oceanbulkavg=.false.
	126	oceanalbvary=.false.
	127	write(,)"Activate runnoff into oceans?"
	128	activerunoff=.false.
[1315]	129	call getin_p("activerunoff",activerunoff)
[875]	130	write(,)" activerunoff = ",activerunoff
	131
[1537]	132	if (activerunoff) then
	133	ALLOCATE(runoff(ngrid))
	134	runoff(1:ngrid)=0
	135	endif
[787]	136
[253]	137	ivap=igcm_h2o_vap
	138	iliq=igcm_h2o_vap
	139	iice=igcm_h2o_ice
	140
	141	write(,) "hydrol: ivap=",ivap
	142	write(,) " iliq=",iliq
	143	write(,) " iice=",iice
	144
	145	! Here's the deal: iice is used in place of igcm_h2o_ice both on the
	146	! surface and in the atmosphere. ivap is used in
	147	! place of igcm_h2o_vap ONLY in the atmosphere, while
	148	! iliq is used in place of igcm_h2o_vap ONLY on the
	149	! surface.
	150	! Soon to be extended to the entire water cycle...
	151
[3100]	152	! LOCAL ocean surface area
[253]	153	oceanarea=0.
[787]	154	do ig=1,ngrid
[1297]	155	if(nint(rnat(ig)).eq.0)then
[1542]	156	oceanarea=oceanarea+cell_area(ig)
[253]	157	endif
	158	enddo
	159
	160	if(oceanbulkavg.and.(oceanarea.le.0.))then
	161	print*,'How are we supposed to average the ocean'
	162	print*,'temperature, when there are no oceans?'
	163	call abort
	164	endif
	165
	166	if(activerunoff.and.(oceanarea.le.0.))then
	167	print*,'You have enabled runoff, but you have no oceans.'
	168	print*,'Where did you think the water was going to go?'
	169	call abort
	170	endif
	171
	172	firstcall = .false.
	173	endif
	174
[3266]	175	! write (,) "oceanarea", oceanarea
[3100]	176
[253]	177	! add physical tendencies already calculated
	178	! ------------------------------------------
	179
[787]	180	do ig=1,ngrid
[253]	181	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
	182	pdtsurf_hyd(ig)=0.0
[787]	183	do iq=1,nq
[253]	184	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
	185	enddo
	186	enddo
	187
[787]	188	do ig=1,ngrid
	189	do iq=1,nq
[253]	190	dqs_hyd(ig,iq) = 0.0
	191	enddo
	192	enddo
	193
[787]	194	do ig = 1, ngrid
[253]	195
[3267]	196	! Ocean regions (rnat = 0)
	197	! -----------------------
[1297]	198	if(nint(rnat(ig)).eq.0)then
[253]	199
[3267]	200	! Parameterization (not used for the moment) to compute the effect of solar zenith angle on the albedo
	201	! --------------------------
	202	!
[1297]	203	! if(diurnal.and.oceanalbvary)then
	204	! fauxo = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
	205	! albedo(ig) = 1.1( .03 + .630/( 1. + fauxofauxo))
	206	! albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
	207	! else
[3267]	208	!
	209	! do nw=1,L_NSPECTV
	210	! albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in inifis_mod.F90. Later we could introduce spectral dependency for alb_ocean.
	211	! enddo
[1297]	212	! end if
[253]	213
[3267]	214	! we first start by fixing the albedo of oceanic grid to that of the ocean
	215	do nw=1,L_NSPECTV
	216	albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in inifis_mod.F90. Later we could introduce spectral dependency for alb_ocean.
	217	enddo
[1297]	218
[3267]	219	if(ok_slab_ocean) then ! if ocean heat transport param activated
[1482]	220
[3268]	221	frac_snow = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/snow_min)) ! Critical snow height (in kg/m2) from ocean_slab_ice routine.
	222	! Standard value should be 15kg/m2 (i.e. about 5 cm). Note that in the previous ocean param. (from BC2014), this value was 45kg/m2 (i.e. about 15cm).
[3291]	223
[1482]	224	! Albedo final calculation :
	225	do nw=1,L_NSPECTV
[3291]	226	alb_ice=albedo_snow_SPECTV(nw)-(albedo_snow_SPECTV(nw)-alb_ice_min)*exp(-sea_ice(ig)/h_alb_ice) ! this replaces the formulation from BC2014
	227	! More details on the parameterization of sea ice albedo vs thickness is provided in the wiki :
	228	! https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Slab_ocean_model
	229	! sea_ice is the ice thickness (calculated in ocean_slab routine) in kg/m2 ; h_alb_ice is fixed to 275.1kg/m2 i.e. 30cm based on comparisons with Brandt et al. 2005
[1482]	230	albedo(ig,nw) = pctsrf_sic(ig)* &
[3268]	231	(albedo_snow_SPECTV(nw)frac_snow + alb_ice(1.0-frac_snow)) &
[1482]	232	+ (1.-pctsrf_sic(ig))*alb_ocean
	233	enddo
[1297]	234
[1482]	235	! Oceanic ice height, just for diagnostics
	236	hice(ig) = MIN(10.,sea_ice(ig)/rhowater)
[3267]	237	else !ok_slab_ocean ; here this is the case where we are dealing with a static ocean
[1297]	238
	239
[3267]	240	! calculate oceanic ice height including the latent heat of ice formation
	241	! hice is the height of oceanic ice with a maximum of maxicethick.
[1482]	242	hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
[3267]	243	twater(ig) = ztsurf(ig) - hice(ig)RLFTTrhowater/pcapcal(ig) ! this is the temperature water would have if we melted the entire ocean ice layer
[1482]	244	hicebis(ig) = hice(ig)
	245	hice(ig) = 0.
[253]	246
[1482]	247	if(twater(ig) .lt. T_h2O_ice_liq)then
	248	E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))pcapcal(ig),RLFTTrhowater*maxicethick)
	249	hice(ig) = E/(RLFTT*rhowater)
	250	hice(ig) = max(hice(ig),0.0)
	251	hice(ig) = min(hice(ig),maxicethick)
	252	pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))RLFTTrhowater/pcapcal(ig)/ptimestep
	253	do nw=1,L_NSPECTV
	254	albedo(ig,nw) = albedo_snow_SPECTV(nw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
	255	enddo
[253]	256
[1482]	257	! if (zqsurf(ig,iice).ge.snowlayer) then
	258	! albedo(ig) = albedoice
	259	! else
	260	! albedo(ig) = albedoocean &
	261	! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
	262	! endif
[253]	263
[1482]	264	else
[253]	265
[1482]	266	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
	267	DO nw=1,L_NSPECTV
	268	albedo(ig,nw) = alb_ocean
	269	ENDDO
[253]	270
[1482]	271	endif
[253]	272
[1482]	273	zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
	274	zqsurf(ig,iice) = hice(ig)*rhowater
[253]	275
[1482]	276	endif!(ok_slab_ocean)
[253]	277
[1297]	278
[3267]	279	! Continental regions (rnat = 1)
	280	! -----------------------
[1297]	281	elseif (nint(rnat(ig)).eq.1) then
[253]	282
[3267]	283	! melt the snow
[650]	284	if(ztsurf(ig).gt.T_h2O_ice_liq)then
[253]	285	if(zqsurf(ig,iice).gt.1.0e-8)then
	286
[650]	287	a = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
[253]	288	b = zqsurf(ig,iice)
	289	fsnoi = min(a,b)
	290
	291	zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
	292	zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi
	293
[3267]	294	! thermal effects
[253]	295	pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep
	296
	297	endif
	298	else
	299
[3267]	300	! freeze the water
[253]	301	if(zqsurf(ig,iliq).gt.1.0e-8)then
	302
[650]	303	a = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
[253]	304	b = zqsurf(ig,iliq)
	305
	306	fsnoi = min(a,b)
	307
	308	zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
	309	zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi
	310
[3267]	311	! thermal effects
[253]	312	pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep
	313
	314	endif
	315	endif
	316
[3267]	317	! add runoff (to simulate transport of water from continental regions to oceanic regions ; in practice, this prevents liquid water to build up too much on continental regions)
[253]	318	if(activerunoff)then
	319	runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
[787]	320	if(ngrid.gt.1)then ! runoff only exists in 3D
[253]	321	if(runoff(ig).ne.0.0)then
	322	zqsurf(ig,iliq) = mx_eau_sol
[3267]	323	! note: runoff is added to ocean at end
[253]	324	endif
	325	end if
	326
	327	endif
	328
[3267]	329	! re-calculate continental albedo
[1482]	330	DO nw=1,L_NSPECTV
	331	albedo(ig,nw) = albedo_bareground(ig)
	332	ENDDO
[253]	333	if (zqsurf(ig,iice).ge.snowlayer) then
[1482]	334	DO nw=1,L_NSPECTV
	335	albedo(ig,nw) = albedo_snow_SPECTV(nw)
	336	ENDDO
[253]	337	else
[1482]	338	DO nw=1,L_NSPECTV
	339	albedo(ig,nw) = albedo_bareground(ig) &
	340	+ (albedo_snow_SPECTV(nw) - albedo_bareground(ig)) &
	341	*zqsurf(ig,iice)/snowlayer
	342	ENDDO
[253]	343	endif
	344
	345	else
	346
	347	print*,'Surface type not recognised in hydrol.F!'
	348	print*,'Exiting...'
	349	call abort
	350
	351	endif
	352
[787]	353	end do ! ig=1,ngrid
[253]	354
[3100]	355
[3266]	356	! simple parameterization to perform crude bulk averaging of temperature in ocean
[253]	357	! ----------------------------------------------------
[3266]	358	if(oceanbulkavg)then
[253]	359
	360	oceanarea2=0.
[787]	361	DO ig=1,ngrid
[1297]	362	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
[1542]	363	oceanarea2=oceanarea2+cell_area(ig)*pcapcal(ig)
[253]	364	end if
	365	END DO
	366
	367	tsea=0.
[787]	368	DO ig=1,ngrid
[1297]	369	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
[1542]	370	tsea=tsea+ztsurf(ig)cell_area(ig)pcapcal(ig)/oceanarea2
[253]	371	end if
	372	END DO
	373
[787]	374	DO ig=1,ngrid
[1297]	375	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0))then
[253]	376	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
	377	end if
	378	END DO
	379
[305]	380	print*,'Mean ocean temperature = ',tsea,' K'
[253]	381
	382	endif
	383
	384	! shove all the runoff water into the ocean
	385	! -----------------------------------------
	386	if(activerunoff)then
	387
[3100]	388	! totalrunoff=0.
[787]	389	do ig=1,ngrid
[3100]	390	runoffamount(ig) = cell_area(ig)*runoff(ig)
	391	! if (nint(rnat(ig)).eq.1) then
	392	! totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
	393	! endif
[253]	394	enddo
[3100]	395
[3267]	396	! collect on the full grid
[3100]	397	call gather(runoffamount,runoffamount_glo)
	398	call gather(zqsurf(1:ngrid,iliq),zqsurf_iliq_glo)
	399	call gather(rnat,rnat_glo)
	400	call gather(cell_area,cell_area_glo)
	401
	402	if (is_master) then
	403	totalrunoff=0.
	404	oceanarea_glo=0.
	405	do ig=1,klon_glo
	406	if (nint(rnat_glo(ig)).eq.1) then
	407	totalrunoff = totalrunoff + runoffamount_glo(ig)
	408	endif
	409	if (nint(rnat_glo(ig)).eq.0) then
	410	oceanarea_glo = oceanarea_glo + cell_area_glo(ig)
	411	endif
	412	enddo
	413
	414	do ig=1,klon_glo
	415	if (nint(rnat_glo(ig)).eq.0) then
	416	zqsurf_iliq_glo(ig) = zqsurf_iliq_glo(ig) + &
	417	totalrunoff/oceanarea_glo
	418	endif
	419	enddo
	420
	421	endif! is_master
	422
[3267]	423	! scatter the field back on all processes
[3100]	424	call scatter(zqsurf_iliq_glo,zqsurf(1:ngrid,iliq))
	425
	426
[253]	427
[3100]	428	! do ig=1,ngrid
	429	! if (nint(rnat(ig)).eq.0) then
	430	! zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
	431	! totalrunoff/oceanarea
	432	! endif
	433	! enddo
[253]	434
[3100]	435	endif !activerunoff
[253]	436
	437	! Re-add the albedo effects of CO2 ice if necessary
	438	! -------------------------------------------------
	439	if(co2cond)then
	440
[787]	441	do ig=1,ngrid
[1482]	442	if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015
	443	DO nw=1,L_NSPECTV
	444	albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
	445	ENDDO
[253]	446	endif
[1482]	447	enddo ! ngrid
[253]	448
[1482]	449	endif ! co2cond
[253]	450
	451
[1484]	452	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 !
	453	dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep - dqsurf(ig,iliq)
	454	dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep - dqsurf(ig,iice)
[253]	455	enddo
	456
[965]	457	if (activerunoff) then
	458	call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
	459	endif
[253]	460
	461	return
	462	end subroutine hydrol

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

Download in other formats:

Original Format