My Project
 All Classes Files Functions Variables Macros
calcul_fluxs_mod.F90
Go to the documentation of this file.
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 
49  include "YOETHF.h"
50  include "FCTTRE.h"
51  include "indicesol.h"
52  include "YOMCST.h"
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 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta
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
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)
250 
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 !
289 END MODULE calcul_fluxs_mod