Context Navigation

calcul_fluxs_mod.F90 @ 1749

Last change on this file since 1749 was 1107, checked in by lguez, 16 years ago
"comconst.h" and "comgeom2.h" are now both fixed and free form. Removed calls to procedure "flush". Corrected kinds of constants which appeared as arguments to "min" or "max" (all arguments are now of the same type and kind).
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 10.9 KB

Rev	Line
[782]	1	!
	2	MODULE calcul_fluxs_mod
	3
	4
	5	CONTAINS
	6	SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
	7	tsurf, p1lay, cal, beta, coef1lay, ps, &
	8	precip_rain, precip_snow, snow, qsurf, &
	9	radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
	10	petAcoef, peqAcoef, petBcoef, peqBcoef, &
	11	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
	12
	13	USE dimphy, ONLY : klon
	14
	15	! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
	16	! une temperature de surface (au cas ou ok_veget = false)
	17	!
	18	! L. Fairhead 4/2000
	19	!
	20	! input:
	21	! knon nombre de points a traiter
	22	! nisurf surface a traiter
	23	! tsurf temperature de surface
	24	! p1lay pression 1er niveau (milieu de couche)
	25	! cal capacite calorifique du sol
	26	! beta evap reelle
	27	! coef1lay coefficient d'echange
	28	! ps pression au sol
	29	! precip_rain precipitations liquides
	30	! precip_snow precipitations solides
	31	! snow champs hauteur de neige
	32	! runoff runoff en cas de trop plein
	33	! petAcoef coeff. A de la resolution de la CL pour t
	34	! peqAcoef coeff. A de la resolution de la CL pour q
	35	! petBcoef coeff. B de la resolution de la CL pour t
	36	! peqBcoef coeff. B de la resolution de la CL pour q
	37	! radsol rayonnement net aus sol (LW + SW)
	38	! dif_grnd coeff. diffusion vers le sol profond
	39	!
	40	! output:
	41	! tsurf_new temperature au sol
	42	! qsurf humidite de l'air au dessus du sol
	43	! fluxsens flux de chaleur sensible
	44	! fluxlat flux de chaleur latente
	45	! dflux_s derivee du flux de chaleur sensible / Ts
	46	! dflux_l derivee du flux de chaleur latente / Ts
	47	!
	48
[793]	49	INCLUDE "YOETHF.h"
	50	INCLUDE "FCTTRE.h"
	51	INCLUDE "indicesol.h"
	52	INCLUDE "YOMCST.h"
[782]	53
	54	! Parametres d'entree
	55	!****************************************************************************************
	56	INTEGER, INTENT(IN) :: knon, nisurf
	57	REAL , INTENT(IN) :: dtime
	58	REAL, DIMENSION(klon), INTENT(IN) :: petAcoef, peqAcoef
	59	REAL, DIMENSION(klon), INTENT(IN) :: petBcoef, peqBcoef
	60	REAL, DIMENSION(klon), INTENT(IN) :: ps, q1lay
	61	REAL, DIMENSION(klon), INTENT(IN) :: tsurf, p1lay, cal, beta, coef1lay
	62	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow ! pas utiles
	63	REAL, DIMENSION(klon), INTENT(IN) :: radsol, dif_grnd
	64	REAL, DIMENSION(klon), INTENT(IN) :: t1lay, u1lay, v1lay
	65
	66	! Parametres entree-sorties
	67	!****************************************************************************************
	68	REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! snow pas utile
	69
	70	! Parametres sorties
	71	!****************************************************************************************
	72	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
	73	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new, evap, fluxsens, fluxlat
	74	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
	75
	76	! Variables locales
	77	!****************************************************************************************
	78	INTEGER :: i
	79	REAL, DIMENSION(klon) :: zx_mh, zx_nh, zx_oh
	80	REAL, DIMENSION(klon) :: zx_mq, zx_nq, zx_oq
	81	REAL, DIMENSION(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
	82	REAL, DIMENSION(klon) :: zx_sl, zx_k1
	83	REAL, DIMENSION(klon) :: d_ts
	84	REAL :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
	85	REAL :: qsat_new, q1_new
	86	REAL, PARAMETER :: t_grnd = 271.35, t_coup = 273.15
	87	REAL, PARAMETER :: max_eau_sol = 150.0
	88	CHARACTER (len = 20) :: modname = 'calcul_fluxs'
	89	LOGICAL :: fonte_neige
	90	LOGICAL, SAVE :: check = .FALSE.
	91	!$OMP THREADPRIVATE(check)
	92
	93	! End definition
	94	!****************************************************************************************
	95
	96	IF (check) WRITE(,)'Entree ', modname,' surface = ',nisurf
	97
	98	IF (check) THEN
	99	WRITE(,)' radsol (min, max)', &
	100	MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
	101	ENDIF
	102
	103	! Traitement neige et humidite du sol
	104	!****************************************************************************************
	105	!
	106	!!$ WRITE(,)'test calcul_flux, surface ', nisurf
	107	!!PB test
	108	!!$ if (nisurf == is_oce) then
	109	!!$ snow = 0.
	110	!!$ qsol = max_eau_sol
	111	!!$ else
	112	!!$ where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
	113	!!$ where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
	114	!!$! snow = max(0.0, snow + (precip_snow - evap) * dtime)
	115	!!$ where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
	116	!!$ endif
	117	!!$ IF (nisurf /= is_ter) qsol = max_eau_sol
	118
	119
	120	!
	121	! Initialisation
	122	!****************************************************************************************
	123	evap = 0.
	124	fluxsens=0.
	125	fluxlat=0.
	126	dflux_s = 0.
	127	dflux_l = 0.
	128	!
	129	! zx_qs = qsat en kg/kg
	130	!****************************************************************************************
	131	DO i = 1, knon
	132	zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
	133	IF (thermcep) THEN
	134	zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
	135	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
	136	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
	137	zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
	138	zx_qs=MIN(0.5,zx_qs)
	139	zcor=1./(1.-retv*zx_qs)
	140	zx_qs=zx_qs*zcor
	141	zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
	142	/RLVTT / zx_pkh(i)
	143	ELSE
	144	IF (tsurf(i).LT.t_coup) THEN
	145	zx_qs = qsats(tsurf(i)) / ps(i)
	146	zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
	147	/ zx_pkh(i)
	148	ELSE
	149	zx_qs = qsatl(tsurf(i)) / ps(i)
	150	zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
	151	/ zx_pkh(i)
	152	ENDIF
	153	ENDIF
	154	zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
	155	zx_qsat(i) = zx_qs
	156	zx_coef(i) = coef1lay(i) * &
	157	(1.0+SQRT(u1lay(i)2+v1lay(i)2)) * &
	158	p1lay(i)/(RD*t1lay(i))
	159
	160	ENDDO
	161
	162
	163	! === Calcul de la temperature de surface ===
	164	! zx_sl = chaleur latente d'evaporation ou de sublimation
	165	!****************************************************************************************
	166
	167	DO i = 1, knon
	168	zx_sl(i) = RLVTT
	169	IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
	170	zx_k1(i) = zx_coef(i)
	171	ENDDO
	172
	173
	174	DO i = 1, knon
	175	! Q
	176	zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
	177	zx_mq(i) = beta(i) * zx_k1(i) * &
	178	(peqAcoef(i) - zx_qsat(i) + &
	179	zx_dq_s_dt(i) * tsurf(i)) &
	180	/ zx_oq(i)
	181	zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
	182	/ zx_oq(i)
	183
	184	! H
	185	zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
	186	zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
	187	zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
	188
	189	! Tsurface
	190	tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
	191	(radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
	192	+ dif_grnd(i) * t_grnd * dtime)/ &
	193	( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
	194	zx_nh(i) + zx_sl(i) * zx_nq(i)) &
	195	+ dtime * dif_grnd(i))
	196
	197	!
	198	! Y'a-t-il fonte de neige?
	199	!
	200	! fonte_neige = (nisurf /= is_oce) .AND. &
	201	! & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
	202	! & .AND. (tsurf_new(i) >= RTT)
	203	! if (fonte_neige) tsurf_new(i) = RTT
	204	d_ts(i) = tsurf_new(i) - tsurf(i)
	205	! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
	206	! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
	207
	208	!== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
	209	!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
	210	evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
	211	fluxlat(i) = - evap(i) * zx_sl(i)
	212	fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
	213
	214	! Derives des flux dF/dTs (W m-2 K-1):
	215	dflux_s(i) = zx_nh(i)
	216	dflux_l(i) = (zx_sl(i) * zx_nq(i))
	217
	218	! Nouvelle valeure de l'humidite au dessus du sol
	219	qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
	220	q1_new = peqAcoef(i) - peqBcoef(i)evap(i)dtime
	221	qsurf(i)=q1_new(1.-beta(i)) + beta(i)qsat_new
	222	!
	223	! en cas de fonte de neige
	224	!
	225	! if (fonte_neige) then
	226	! bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
	227	! & dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
	228	! & RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
	229	! bilan_f = max(0., bilan_f)
	230	! fq_fonte = bilan_f / zx_sl(i)
	231	! snow(i) = max(0., snow(i) - fq_fonte * dtime)
	232	! qsol(i) = qsol(i) + (fq_fonte * dtime)
	233	! endif
	234	!!$ if (nisurf == is_ter) &
	235	!!$ & run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
	236	!!$ qsol(i) = min(qsol(i), max_eau_sol)
	237	ENDDO
	238	!
	239	!****************************************************************************************
	240	!
	241	END SUBROUTINE calcul_fluxs
[1067]	242	!
	243	!****************************************************************************************
	244	!
	245	SUBROUTINE calcul_flux_wind(knon, dtime, &
	246	u0, v0, u1, v1, cdrag_m, &
	247	AcoefU, AcoefV, BcoefU, BcoefV, &
	248	p1lay, t1lay, &
	249	flux_u1, flux_v1)
[782]	250
[1067]	251	USE dimphy
	252	INCLUDE "YOMCST.h"
	253
	254	! Input arguments
	255	!****************************************************************************************
	256	INTEGER, INTENT(IN) :: knon
	257	REAL, INTENT(IN) :: dtime
	258	REAL, DIMENSION(klon), INTENT(IN) :: u0, v0 ! u and v at niveau 0
	259	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1 ! u and v at niveau 1
	260	REAL, DIMENSION(klon), INTENT(IN) :: cdrag_m ! cdrag pour momentum
	261	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
	262	REAL, DIMENSION(klon), INTENT(IN) :: p1lay ! pression 1er niveau (milieu de couche)
	263	REAL, DIMENSION(klon), INTENT(IN) :: t1lay ! temperature
	264	! Output arguments
	265	!****************************************************************************************
	266	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1
	267	REAL, DIMENSION(klon), INTENT(OUT) :: flux_v1
	268
	269	! Local variables
	270	!****************************************************************************************
	271	INTEGER :: i
	272	REAL :: mod_wind, buf
	273
	274	!****************************************************************************************
	275	! Calculate the surface flux
	276	!
	277	!****************************************************************************************
	278	DO i=1,knon
	279	mod_wind = 1.0 + SQRT((u1(i) - u0(i))2 + (v1(i)-v0(i))2)
	280	buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
	281	flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
	282	flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
	283	END DO
	284
	285	END SUBROUTINE calcul_flux_wind
	286	!
	287	!****************************************************************************************
	288	!
[782]	289	END MODULE calcul_fluxs_mod

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

Download in other formats:

Original Format