Context Navigation

surf_landice_mod.F90 @ 3851

Last change on this file since 3851 was 3851, checked in by dcugnet, 4 years ago
Update the branch to the current trunk.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 18.0 KB

Rev	Line
[781]	1	!
	2	MODULE surf_landice_mod
	3
	4	IMPLICIT NONE
	5
	6	CONTAINS
	7	!
	8	!****************************************************************************************
	9	!
	10	SUBROUTINE surf_landice(itime, dtime, knon, knindex, &
[1865]	11	rlon, rlat, debut, lafin, &
	12	rmu0, lwdownm, albedo, pphi1, &
[888]	13	swnet, lwnet, tsurf, p1lay, &
[1067]	14	cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
	15	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	16	AcoefU, AcoefV, BcoefU, BcoefV, &
[2240]	17	ps, u1, v1, gustiness, rugoro, pctsrf, &
[888]	18	snow, qsurf, qsol, agesno, &
[2243]	19	tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, &
[1067]	20	tsurf_new, dflux_s, dflux_l, &
[1865]	21	slope, cloudf, &
	22	snowhgt, qsnow, to_ice, sissnow, &
	23	alb3, runoff, &
[1067]	24	flux_u1, flux_v1)
[781]	25
[1067]	26	USE dimphy
[3851]	27	USE surface_data, ONLY : type_ocean, calice, calsno, landice_opt, n_dtis
[1067]	28	USE fonte_neige_mod, ONLY : fonte_neige, run_off_lic
	29	USE cpl_mod, ONLY : cpl_send_landice_fields
	30	USE calcul_fluxs_mod
[1334]	31	USE phys_output_var_mod
[2728]	32	!FC
	33	USE ioipsl_getin_p_mod, ONLY : getin_p
	34
[1865]	35	#ifdef CPP_SISVAT
	36	USE surf_sisvat_mod, ONLY : surf_sisvat
	37	#endif
[3851]	38
	39	#ifdef CPP_INLANDSIS
	40	USE surf_inlandsis_mod, ONLY : surf_inlandsis
	41	#endif
	42
[1785]	43	USE indice_sol_mod
[1067]	44
[1785]	45	! INCLUDE "indicesol.h"
[781]	46	INCLUDE "dimsoil.h"
[793]	47	INCLUDE "YOMCST.h"
	48	INCLUDE "clesphys.h"
[781]	49
	50	! Input variables
	51	!****************************************************************************************
	52	INTEGER, INTENT(IN) :: itime, knon
	53	INTEGER, DIMENSION(klon), INTENT(in) :: knindex
	54	REAL, INTENT(in) :: dtime
[888]	55	REAL, DIMENSION(klon), INTENT(IN) :: swnet ! net shortwave radiance
	56	REAL, DIMENSION(klon), INTENT(IN) :: lwnet ! net longwave radiance
[781]	57	REAL, DIMENSION(klon), INTENT(IN) :: tsurf
	58	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
[1067]	59	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
[781]	60	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
	61	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
[1067]	62	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
	63	REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
	64	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
[781]	65	REAL, DIMENSION(klon), INTENT(IN) :: ps
[2240]	66	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
[781]	67	REAL, DIMENSION(klon), INTENT(IN) :: rugoro
	68	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
	69
[1865]	70	LOGICAL, INTENT(IN) :: debut !true if first step
	71	LOGICAL, INTENT(IN) :: lafin !true if last step
	72	REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
	73	REAL, DIMENSION(klon), INTENT(IN) :: rmu0
	74	REAL, DIMENSION(klon), INTENT(IN) :: lwdownm !ylwdown
	75	REAL, DIMENSION(klon), INTENT(IN) :: albedo !mean albedo
	76	REAL, DIMENSION(klon), INTENT(IN) :: pphi1
	77	REAL, DIMENSION(klon), INTENT(IN) :: slope !mean slope in grid box
	78	REAL, DIMENSION(klon), INTENT(IN) :: cloudf !total cloud fraction
	79
[781]	80	! In/Output variables
	81	!****************************************************************************************
	82	REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
	83	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
	84	REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
	85
	86	! Output variables
	87	!****************************************************************************************
	88	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
[2243]	89	REAL, DIMENSION(klon), INTENT(OUT) :: z0m, z0h
[2227]	90	!albedo SB >>>
	91	! REAL, DIMENSION(klon), INTENT(OUT) :: alb1 ! new albedo in visible SW interval
	92	! REAL, DIMENSION(klon), INTENT(OUT) :: alb2 ! new albedo in near IR interval
[3851]	93	REAL, DIMENSION(6), INTENT(IN) :: SFRWL
	94	REAL, DIMENSION(klon,nsw), INTENT(OUT) :: alb_dir,alb_dif
[2227]	95	!albedo SB <<<
[781]	96	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
[888]	97	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
[781]	98	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
[1067]	99	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
[781]	100
[1865]	101	REAL, DIMENSION(klon), INTENT(OUT) :: alb3
	102	REAL, DIMENSION(klon), INTENT(OUT) :: qsnow !column water in snow [kg/m2]
	103	REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt !Snow height (m)
	104	REAL, DIMENSION(klon), INTENT(OUT) :: to_ice
	105	REAL, DIMENSION(klon), INTENT(OUT) :: sissnow
	106	REAL, DIMENSION(klon), INTENT(OUT) :: runoff !Land ice runoff
	107
	108
[781]	109	! Local variables
	110	!****************************************************************************************
	111	REAL, DIMENSION(klon) :: soilcap, soilflux
	112	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
	113	REAL, DIMENSION(klon) :: zfra, alb_neig
[888]	114	REAL, DIMENSION(klon) :: radsol
[3851]	115	REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay, ustar
	116	INTEGER :: i,j,nt
[781]	117
[1865]	118	REAL, DIMENSION(klon) :: emis_new !Emissivity
	119	REAL, DIMENSION(klon) :: swdown,lwdown
	120	REAL, DIMENSION(klon) :: precip_snow_adv, snow_adv !Snow Drift precip./advection
[3851]	121	REAL, DIMENSION(klon) :: zsl_height, wind_velo !surface layer height, wind spd
[1865]	122	REAL, DIMENSION(klon) :: dens_air, snow_cont_air !air density; snow content air
	123	REAL, DIMENSION(klon) :: alb_soil !albedo of underlying ice
	124	REAL, DIMENSION(klon) :: pexner !Exner potential
	125	REAL :: pref
[3851]	126	REAL, DIMENSION(klon,nsoilmx) :: tsoil0 !modif
	127	REAL :: dtis ! subtimestep
	128	LOGICAL :: debut_is, lafin_is ! debut and lafin for inlandsis
[1865]	129
	130	CHARACTER (len = 20) :: modname = 'surf_landice'
	131	CHARACTER (len = 80) :: abort_message
	132
[2728]	133
[2227]	134	!albedo SB >>>
	135	real,dimension(klon) :: alb1,alb2
	136	!albedo SB <<<
	137
[781]	138	! End definition
	139	!****************************************************************************************
[2728]	140	!FC
	141	!FC
	142	REAL,SAVE :: alb_vis_sno_lic
	143	!$OMP THREADPRIVATE(alb_vis_sno_lic)
	144	REAL,SAVE :: alb_nir_sno_lic
	145	!$OMP THREADPRIVATE(alb_nir_sno_lic)
	146	LOGICAL, SAVE :: firstcall = .TRUE.
	147	!$OMP THREADPRIVATE(firstcall)
	148
	149
[3851]	150	!FC firtscall initializations
	151	!******************************************************************************************
[2728]	152	IF (firstcall) THEN
	153	alb_vis_sno_lic=0.77
	154	CALL getin_p('alb_vis_sno_lic',alb_vis_sno_lic)
	155	PRINT*, 'alb_vis_sno_lic',alb_vis_sno_lic
	156	alb_nir_sno_lic=0.77
	157	CALL getin_p('alb_nir_sno_lic',alb_nir_sno_lic)
	158	PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
[3851]	159
	160	! z0m=1.e-3
	161	! z0h = z0m
[2728]	162	firstcall=.false.
	163	ENDIF
[3851]	164	!******************************************************************************************
	165
[781]	166	! Initialize output variables
[1865]	167	alb3(:) = 999999.
[888]	168	alb2(:) = 999999.
	169	alb1(:) = 999999.
[1865]	170
	171	runoff(:) = 0.
[888]	172	!****************************************************************************************
	173	! Calculate total absorbed radiance at surface
	174	!
	175	!****************************************************************************************
	176	radsol(:) = 0.0
	177	radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
[781]	178
	179	!****************************************************************************************
[3851]	180	! landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ...
	181	! landice_opt = 1 : prepare and call SISVAT snow model
	182	! landice_opt = 2 : prepare and call INLANDSIS snow model
[1865]	183	!****************************************************************************************
[3851]	184
	185
	186	IF (landice_opt .EQ. 1) THEN
	187
	188	!****************************************************************************************
	189	! CALL to SISVAT interface
	190	!****************************************************************************************
	191
[1865]	192	#ifdef CPP_SISVAT
	193	! Prepare for calling SISVAT
	194
	195	! Calculate incoming flux for SW and LW interval: swdown, lwdown
	196	swdown(:) = 0.0
	197	lwdown(:) = 0.0
	198	DO i = 1, knon
	199	swdown(i) = swnet(i)/(1-albedo(i))
	200	lwdown(i) = lwdownm(i)
	201	END DO
	202
	203	! Set constants and compute some input for SISVAT
	204	snow_adv(:) = 0. ! no snow blown in for now
	205	snow_cont_air(:) = 0.
	206	alb_soil(:) = albedo(:)
	207	pref = 100000. ! = 1000 hPa
	208	DO i = 1, knon
	209	wind_velo(i) = u1(i)2 + v1(i)2
	210	wind_velo(i) = wind_velo(i)**0.5
	211	pexner(i) = (p1lay(i)/pref)**(RD/RCPD)
	212	dens_air(i) = p1lay(i)/RD/temp_air(i) ! dry air density
[3851]	213	zsl_height(i) = pphi1(i)/RG
[1865]	214	END DO
	215
[3851]	216
[1865]	217	! config: compute everything with SV but temperatures afterwards with soil/calculfluxs
	218	DO i = 1, knon
	219	tsoil0(i,:)=tsoil(i,:)
	220	END DO
[2254]	221	! Martin
	222	PRINT*, 'on appelle surf_sisvat'
	223	! Martin
[1865]	224	CALL surf_sisvat(knon, rlon, rlat, knindex, itime, dtime, debut, lafin, &
	225	rmu0, swdown, lwdown, pexner, ps, p1lay, &
	226	precip_rain, precip_snow, precip_snow_adv, snow_adv, &
[3851]	227	zsl_height, wind_velo, temp_air, dens_air, spechum, tsurf, &
[1865]	228	rugoro, snow_cont_air, alb_soil, slope, cloudf, &
	229	radsol, qsol, tsoil0, snow, snowhgt, qsnow, to_ice,sissnow, agesno, &
	230	AcoefH, AcoefQ, BcoefH, BcoefQ, cdragh, &
	231	run_off_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &
	232	tsurf_new, alb1, alb2, alb3, &
[2243]	233	emis_new, z0m, qsurf)
	234	z0h(1:knon)=z0m(1:knon) ! en attendant mieux
[1865]	235
	236	! Suppose zero surface speed
	237	u0(:) = 0.0
	238	v0(:) = 0.0
	239	! The calculation of heat/water fluxes, otherwise done by "CALL calcul_fluxs" is
	240	! integrated in SISVAT, using the same method. It can be found in "sisvat.f", in the
	241	! subroutine "SISVAT_TS2".
	242	! u0, v0=0., dif_grnd=0. and beta=1 are assumed there!
	243
	244	CALL calcul_flux_wind(knon, dtime, &
[2240]	245	u0, v0, u1, v1, gustiness, cdragm, &
[1865]	246	AcoefU, AcoefV, BcoefU, BcoefV, &
	247	p1lay, temp_air, &
	248	flux_u1, flux_v1)
	249	#else
[3851]	250	abort_message='Pb de coherence: landice_opt = 1 mais CPP_SISVAT = .false.'
[2311]	251	CALL abort_physic(modname,abort_message,1)
[1865]	252	#endif
	253
	254	!****************************************************************************************
[3851]	255	! CALL to INLANDSIS interface
	256	!****************************************************************************************
	257
	258	ELSE IF (landice_opt .EQ. 2) THEN
	259	#ifdef CPP_INLANDSIS
	260
	261	debut_is=debut
	262	lafin_is=.false.
	263	! Suppose zero surface speed
	264	u0(:) = 0.0
	265	v0(:) = 0.0
	266
	267
	268	CALL calcul_flux_wind(knon, dtime, &
	269	u0, v0, u1, v1, gustiness, cdragm, &
	270	AcoefU, AcoefV, BcoefU, BcoefV, &
	271	p1lay, temp_air, &
	272	flux_u1, flux_v1)
	273
	274
	275	! Set constants and compute some input for SISVAT
	276	! = 1000 hPa
	277	! and calculate incoming flux for SW and LW interval: swdown, lwdown
	278	swdown(:) = 0.0
	279	lwdown(:) = 0.0
	280	snow_adv(:) = 0. ! no snow blown in for now
	281	snow_cont_air(:) = 0.
	282	alb_soil(:) = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
	283	ustar(:) = 0.
	284	pref = 100000.
	285	DO i = 1, knon
	286	swdown(i) = swnet(i)/(1-albedo(i))
	287	lwdown(i) = lwdownm(i)
	288	wind_velo(i) = u1(i)2 + v1(i)2
	289	wind_velo(i) = wind_velo(i)**0.5
	290	pexner(i) = (p1lay(i)/pref)**(RD/RCPD)
	291	dens_air(i) = p1lay(i)/RD/temp_air(i) ! dry air density
	292	zsl_height(i) = pphi1(i)/RG
	293	tsoil0(i,:) = tsoil(i,:)
	294	ustar(i)= (cdragm(i)(wind_velo(i)2))*0.5
	295	END DO
	296
	297
	298	! Subtimestepping
	299
	300	dtis=dtime/n_dtis
	301
	302	DO nt=1,n_dtis
	303
	304	IF (lafin .and. nt.eq.n_dtis) THEN
	305	lafin_is=.true.
	306	END IF
	307
	308	!PRINT*,'RENTRE DANS INLANDSIS','itime',itime,'dtime',dtime,'dtis',dtis
	309	CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is, &
	310	rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, &
	311	precip_rain, precip_snow, precip_snow_adv, snow_adv, &
	312	zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
	313	rugoro, snow_cont_air, alb_soil, slope, cloudf, &
	314	radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice,sissnow, agesno, &
	315	AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
	316	run_off_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &
	317	tsurf_new, alb1, alb2, alb3, &
	318	emis_new, z0m, z0h, qsurf)
	319
	320	debut_is=.false.
	321
	322	END DO
	323
	324
	325	#else
	326	abort_message='Pb de coherence: landice_opt = 2 mais CPP_INLANDSIS = .false.'
	327	CALL abort_physic(modname,abort_message,1)
	328	#endif
	329
	330
	331
	332	ELSE
	333
	334	!****************************************************************************************
[781]	335	! Soil calculations
	336	!
	337	!****************************************************************************************
[3851]	338
	339	! EV: use calbeta
	340	CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
	341
	342
	343	! use soil model and recalculate properly cal
[781]	344	IF (soil_model) THEN
	345	CALL soil(dtime, is_lic, knon, snow, tsurf, tsoil, soilcap, soilflux)
	346	cal(1:knon) = RCPD / soilcap(1:knon)
	347	radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
	348	ELSE
	349	cal = RCPD * calice
	350	WHERE (snow > 0.0) cal = RCPD * calsno
	351	ENDIF
	352
	353
	354	!****************************************************************************************
	355	! Calulate fluxes
	356	!
	357	!****************************************************************************************
[3851]	358	! beta(:) = 1.0
	359	! dif_grnd(:) = 0.0
[781]	360
[1067]	361	! Suppose zero surface speed
	362	u0(:)=0.0
	363	v0(:)=0.0
	364	u1_lay(:) = u1(:) - u0(:)
	365	v1_lay(:) = v1(:) - v0(:)
	366
[781]	367	CALL calcul_fluxs(knon, is_lic, dtime, &
[2254]	368	tsurf, p1lay, cal, beta, cdragh, cdragh, ps, &
[781]	369	precip_rain, precip_snow, snow, qsurf, &
[2240]	370	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
[2254]	371	1.,AcoefH, AcoefQ, BcoefH, BcoefQ, &
[781]	372	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
	373
[1067]	374	CALL calcul_flux_wind(knon, dtime, &
[2240]	375	u0, v0, u1, v1, gustiness, cdragm, &
[1067]	376	AcoefU, AcoefV, BcoefU, BcoefV, &
	377	p1lay, temp_air, &
	378	flux_u1, flux_v1)
[781]	379
	380	!****************************************************************************************
	381	! Calculate snow height, age, run-off,..
	382	!
	383	!****************************************************************************************
[3851]	384	CALL fonte_neige(knon, is_lic, knindex, dtime, &
[781]	385	tsurf, precip_rain, precip_snow, &
	386	snow, qsol, tsurf_new, evap)
	387
	388
	389	!****************************************************************************************
	390	! Calculate albedo
	391	!
	392	!****************************************************************************************
	393	CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
[3851]	394
[781]	395	WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
	396	zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
[888]	397	alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
[781]	398	0.6 * (1.0-zfra(1:knon))
	399	!
	400	!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
[888]	401	! alb1(1 : knon) = 0.6 !IM cf FH/GK
	402	! alb1(1 : knon) = 0.82
	403	! alb1(1 : knon) = 0.77 !211003 Ksta0.77
	404	! alb1(1 : knon) = 0.8 !KstaTER0.8 & LMD_ARMIP5
[781]	405	!IM: KstaTER0.77 & LMD_ARMIP6
	406
[3851]	407	! Attantion: alb1 and alb2 are not the same!
[2728]	408	alb1(1:knon) = alb_vis_sno_lic
	409	alb2(1:knon) = alb_nir_sno_lic
[781]	410
	411
	412	!****************************************************************************************
	413	! Rugosity
	414	!
	415	!****************************************************************************************
[2243]	416	z0m=1.e-3
	417	z0h = z0m
	418	z0m = SQRT(z0m2+rugoro2)
	419
[781]	420
[1865]	421
[3851]	422
	423
	424
	425	END IF ! landice_opt
	426
	427
[781]	428	!****************************************************************************************
	429	! Send run-off on land-ice to coupler if coupled ocean.
[1865]	430	! run_off_lic has been calculated in fonte_neige or surf_sisvat
[781]	431	!
	432	!****************************************************************************************
[996]	433	IF (type_ocean=='couple') THEN
[781]	434	CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic)
	435	ENDIF
[1865]	436
	437	! transfer runoff rate [kg/m2/s](!) to physiq for output
	438	runoff(1:knon)=run_off_lic(1:knon)/dtime
	439
[1334]	440
	441	!****************************************************************************************
[3851]	442	! Etienne: comment these lines because of duplication just below
	443	! snow_o=0.
	444	! zfra_o = 0.
	445	! DO j = 1, knon
	446	! i = knindex(j)
	447	! snow_o(i) = snow(j)
	448	! zfra_o(i) = zfra(j)
	449	! ENDDO
	450	!
[1403]	451	!****************************************************************************************
	452	snow_o=0.
	453	zfra_o = 0.
	454	DO j = 1, knon
	455	i = knindex(j)
	456	snow_o(i) = snow(j)
	457	zfra_o(i) = zfra(j)
	458	ENDDO
	459
	460
[2227]	461	!albedo SB >>>
	462	select case(NSW)
	463	case(2)
	464	alb_dir(1:knon,1)=alb1(1:knon)
	465	alb_dir(1:knon,2)=alb2(1:knon)
	466	case(4)
	467	alb_dir(1:knon,1)=alb1(1:knon)
	468	alb_dir(1:knon,2)=alb2(1:knon)
	469	alb_dir(1:knon,3)=alb2(1:knon)
	470	alb_dir(1:knon,4)=alb2(1:knon)
	471	case(6)
	472	alb_dir(1:knon,1)=alb1(1:knon)
	473	alb_dir(1:knon,2)=alb1(1:knon)
	474	alb_dir(1:knon,3)=alb1(1:knon)
	475	alb_dir(1:knon,4)=alb2(1:knon)
	476	alb_dir(1:knon,5)=alb2(1:knon)
	477	alb_dir(1:knon,6)=alb2(1:knon)
	478	end select
	479	alb_dif=alb_dir
	480	!albedo SB <<<
	481
	482
	483
	484
[781]	485	END SUBROUTINE surf_landice
	486	!
	487	!****************************************************************************************
	488	!
	489	END MODULE surf_landice_mod
	490
	491
	492

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

Download in other formats:

Original Format