source: LMDZ5/branches/testing/libf/phylmd/calcul_fluxs_mod.F90 @ 5403

Last change on this file since 5403 was 2542, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2487:2541 into testing branch
  • 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:executable set to *
  • Property svn:keywords set to Id
File size: 13.7 KB
Line 
1!
2! $Id: calcul_fluxs_mod.F90 2542 2016-06-06 14:04:57Z abarral $
3!
4MODULE calcul_fluxs_mod
5
6
7CONTAINS
8  SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
9       tsurf, p1lay, cal, beta, cdragh, cdragq, ps, &
10       precip_rain, precip_snow, snow, qsurf, &
11       radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, &
12       fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, &
13       tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
14       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
15 
16   
17    USE dimphy, ONLY : klon
18    USE indice_sol_mod
19
20    INCLUDE "clesphys.h"
21
22! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
23! une temperature de surface (au cas ou ok_veget = false)
24!
25! L. Fairhead 4/2000
26!
27! input:
28!   knon         nombre de points a traiter
29!   nisurf       surface a traiter
30!   tsurf        temperature de surface
31!   p1lay        pression 1er niveau (milieu de couche)
32!   cal          capacite calorifique du sol
33!   beta         evap reelle
34!   cdragh       coefficient d'echange temperature
35!   cdragq       coefficient d'echange evaporation
36!   ps           pression au sol
37!   precip_rain  precipitations liquides
38!   precip_snow  precipitations solides
39!   snow         champs hauteur de neige
40!   runoff       runoff en cas de trop plein
41!   petAcoef     coeff. A de la resolution de la CL pour t
42!   peqAcoef     coeff. A de la resolution de la CL pour q
43!   petBcoef     coeff. B de la resolution de la CL pour t
44!   peqBcoef     coeff. B de la resolution de la CL pour q
45!   radsol       rayonnement net aus sol (LW + SW)
46!   dif_grnd     coeff. diffusion vers le sol profond
47!
48! output:
49!   tsurf_new    temperature au sol
50!   qsurf        humidite de l'air au dessus du sol
51!   fluxsens     flux de chaleur sensible
52!   fluxlat      flux de chaleur latente
53!   dflux_s      derivee du flux de chaleur sensible / Ts
54!   dflux_l      derivee du flux de chaleur latente  / Ts
55!   sens_prec_liq flux sensible lié aux echanges de precipitations liquides
56!   sens_prec_sol                                    precipitations solides
57!   lat_prec_liq  flux latent lié aux echanges de precipitations liquides
58!   lat_prec_sol                                  precipitations solides
59
60    INCLUDE "YOETHF.h"
61    INCLUDE "FCTTRE.h"
62    INCLUDE "YOMCST.h"
63
64! Parametres d'entree
65!****************************************************************************************
66    INTEGER, INTENT(IN)                  :: knon, nisurf
67    REAL   , INTENT(IN)                  :: dtime
68    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
69    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
70    REAL, DIMENSION(klon), INTENT(IN)    :: ps, q1lay
71    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf, p1lay, cal, beta, cdragh,cdragq
72    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow ! pas utiles
73    REAL, DIMENSION(klon), INTENT(IN)    :: radsol, dif_grnd
74    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay, u1lay, v1lay,gustiness
75    REAL,                  INTENT(IN)    :: fqsat ! correction factor on qsat (generally 0.98 over salty water, 1 everywhere else)
76
77! Parametres entree-sorties
78!****************************************************************************************
79    REAL, DIMENSION(klon), INTENT(INOUT) :: snow  ! snow pas utile
80
81! Parametres sorties
82!****************************************************************************************
83    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
84    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new, evap, fluxsens, fluxlat
85    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l
86    REAL, DIMENSION(klon), OPTIONAL      :: sens_prec_liq, sens_prec_sol
87    REAL, DIMENSION(klon), OPTIONAL      :: lat_prec_liq, lat_prec_sol
88
89! Variables locales
90!****************************************************************************************
91    INTEGER                              :: i
92    REAL, DIMENSION(klon)                :: zx_mh, zx_nh, zx_oh
93    REAL, DIMENSION(klon)                :: zx_mq, zx_nq, zx_oq
94    REAL, DIMENSION(klon)                :: zx_pkh, zx_dq_s_dt, zx_qsat
95    REAL, DIMENSION(klon)                :: zx_sl, zx_coefh, zx_coefq, zx_wind
96    REAL, DIMENSION(klon)                :: d_ts
97    REAL                                 :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
98    REAL                                 :: qsat_new, q1_new
99    REAL, PARAMETER                      :: t_grnd = 271.35, t_coup = 273.15
100    REAL, PARAMETER                      :: max_eau_sol = 150.0
101    CHARACTER (len = 20)                 :: modname = 'calcul_fluxs'
102    LOGICAL                              :: fonte_neige
103    LOGICAL, SAVE                        :: check = .FALSE.
104    !$OMP THREADPRIVATE(check)
105
106! End definition
107!****************************************************************************************
108
109    IF (check) WRITE(*,*)'Entree ', modname,' surface = ',nisurf
110   
111    IF (check) THEN
112       WRITE(*,*)' radsol (min, max)', &
113            MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
114    ENDIF
115 
116! Traitement neige et humidite du sol
117!****************************************************************************************
118!
119!!$  WRITE(*,*)'test calcul_flux, surface ', nisurf
120!!PB test
121!!$    if (nisurf == is_oce) then
122!!$      snow = 0.
123!!$      qsol = max_eau_sol
124!!$    else
125!!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
126!!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
127!!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)
128!!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
129!!$    endif
130!!$    IF (nisurf /= is_ter) qsol = max_eau_sol
131
132
133!
134! Initialisation
135!****************************************************************************************
136    evap = 0.
137    fluxsens=0.
138    fluxlat=0.
139    dflux_s = 0.
140    dflux_l = 0.
141    if (PRESENT(sens_prec_liq)) sens_prec_liq = 0.
142    if (PRESENT(sens_prec_sol)) sens_prec_sol = 0.
143    if (PRESENT(lat_prec_liq)) lat_prec_liq = 0.
144    if (PRESENT(lat_prec_sol)) lat_prec_sol = 0.
145!
146! zx_qs = qsat en kg/kg
147!****************************************************************************************
148    DO i = 1, knon
149       zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
150       IF (thermcep) THEN
151          zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
152          zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
153          zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
154          zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
155          zx_qs=MIN(0.5,zx_qs)
156          zcor=1./(1.-retv*zx_qs)
157          zx_qs=zx_qs*zcor
158          zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
159               /RLVTT / zx_pkh(i)
160       ELSE
161          IF (tsurf(i).LT.t_coup) THEN
162             zx_qs = qsats(tsurf(i)) / ps(i)
163             zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
164                  / zx_pkh(i)
165          ELSE
166             zx_qs = qsatl(tsurf(i)) / ps(i)
167             zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
168                  / zx_pkh(i)
169          ENDIF
170       ENDIF
171       zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
172       zx_qsat(i) = zx_qs
173       zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)**2+v1lay(i)**2)
174       zx_coefh(i) = cdragh(i) * zx_wind(i) * p1lay(i)/(RD*t1lay(i))
175       zx_coefq(i) = cdragq(i) * zx_wind(i) * p1lay(i)/(RD*t1lay(i))
176!      zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)**2+v1lay(i)**2) &
177!                * p1lay(i)/(RD*t1lay(i))
178!      zx_coefh(i) = cdragh(i) * zx_wind(i)
179!      zx_coefq(i) = cdragq(i) * zx_wind(i)
180    ENDDO
181
182
183! === Calcul de la temperature de surface ===
184! zx_sl = chaleur latente d'evaporation ou de sublimation
185!****************************************************************************************
186
187    DO i = 1, knon
188       zx_sl(i) = RLVTT
189       IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
190    ENDDO
191   
192
193    DO i = 1, knon
194! Q
195       zx_oq(i) = 1. - (beta(i) * zx_coefq(i) * peqBcoef(i) * dtime)
196       zx_mq(i) = beta(i) * zx_coefq(i) * &
197            (peqAcoef(i) -             &
198! conv num avec precedente version
199            fqsat * zx_qsat(i) + fqsat * zx_dq_s_dt(i) * tsurf(i))  &
200!           fqsat * ( zx_qsat(i) - zx_dq_s_dt(i) * tsurf(i)) ) &
201            / zx_oq(i)
202       zx_nq(i) = beta(i) * zx_coefq(i) * (- fqsat * zx_dq_s_dt(i)) &
203            / zx_oq(i)
204       
205! H
206       zx_oh(i) = 1. - (zx_coefh(i) * petBcoef(i) * dtime)
207       zx_mh(i) = zx_coefh(i) * petAcoef(i) / zx_oh(i)
208       zx_nh(i) = - (zx_coefh(i) * RCPD * zx_pkh(i))/ zx_oh(i)
209     
210! Tsurface
211       tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
212            (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
213            + dif_grnd(i) * t_grnd * dtime)/ &
214            ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
215            zx_nh(i) + zx_sl(i) * zx_nq(i)) & 
216            + dtime * dif_grnd(i))
217
218!
219! Y'a-t-il fonte de neige?
220!
221!    fonte_neige = (nisurf /= is_oce) .AND. &
222!     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
223!     & .AND. (tsurf_new(i) >= RTT)
224!    if (fonte_neige) tsurf_new(i) = RTT 
225       d_ts(i) = tsurf_new(i) - tsurf(i)
226!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
227!    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
228
229!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
230!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
231       evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
232       fluxlat(i) = - evap(i) * zx_sl(i)
233       fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
234       
235! Derives des flux dF/dTs (W m-2 K-1):
236       dflux_s(i) = zx_nh(i)
237       dflux_l(i) = (zx_sl(i) * zx_nq(i))
238
239! Nouvelle valeure de l'humidite au dessus du sol
240       qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
241       q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
242       qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
243!
244! en cas de fonte de neige
245!
246!    if (fonte_neige) then
247!      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
248!     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
249!     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
250!      bilan_f = max(0., bilan_f)
251!      fq_fonte = bilan_f / zx_sl(i)
252!      snow(i) = max(0., snow(i) - fq_fonte * dtime)
253!      qsol(i) = qsol(i) + (fq_fonte * dtime)
254!    endif
255!!$    if (nisurf == is_ter)  &
256!!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
257!!$    qsol(i) = min(qsol(i), max_eau_sol)
258!
259! calcul de l'enthalpie des precipitations liquides et solides
260!
261!       if (PRESENT(enth_prec_liq))                   &
262!       enth_prec_liq(i) = rcw * (t1lay(i) - tsurf(i)) * &
263!                          precip_rain(i)
264!       if (PRESENT(enth_prec_sol))                  &
265!       enth_prec_sol(i) = rcs * (t1lay(i) - tsurf(i)) * &
266!                          precip_snow(i)
267! On calcule par rapport a T=0
268       if (PRESENT(sens_prec_liq))                   &
269         sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * &
270                          precip_rain(i)
271       if (PRESENT(sens_prec_sol))                   &
272         sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * &
273                          precip_snow(i)
274       if (PRESENT(lat_prec_liq))                    &
275         lat_prec_liq(i) =  precip_rain(i) * (RLVTT - RLVTT)
276       if (PRESENT(lat_prec_sol))                    &
277         lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
278    ENDDO
279
280   
281!       if (PRESENT(sens_prec_liq))                  &
282!          WRITE(*,*)' calculs_fluxs sens_prec_liq (min, max)', &
283!          MINVAL(sens_prec_liq(1:knon)), MAXVAL(sens_prec_liq(1:knon))
284!       if (PRESENT(sens_prec_sol))                  &
285!          WRITE(*,*)' calculs_fluxs sens_prec_sol (min, max)', &
286!          MINVAL(sens_prec_sol(1:knon)), MAXVAL(sens_prec_sol(1:knon))
287 
288!
289!****************************************************************************************
290!
291  END SUBROUTINE calcul_fluxs
292!
293!****************************************************************************************
294!
295  SUBROUTINE calcul_flux_wind(knon, dtime, &
296       u0, v0, u1, v1, gustiness, cdrag_m, &
297       AcoefU, AcoefV, BcoefU, BcoefV, &
298       p1lay, t1lay, &
299       flux_u1, flux_v1)
300
301    USE dimphy
302    INCLUDE "YOMCST.h"
303    INCLUDE "clesphys.h"
304
305! Input arguments
306!****************************************************************************************
307    INTEGER, INTENT(IN)                  :: knon
308    REAL, INTENT(IN)                     :: dtime
309    REAL, DIMENSION(klon), INTENT(IN)    :: u0, v0  ! u and v at niveau 0
310    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness  ! u and v at niveau 1
311    REAL, DIMENSION(klon), INTENT(IN)    :: cdrag_m ! cdrag pour momentum
312    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
313    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay   ! pression 1er niveau (milieu de couche)
314    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay   ! temperature
315! Output arguments
316!****************************************************************************************
317    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1
318    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v1
319
320! Local variables
321!****************************************************************************************
322    INTEGER                              :: i
323    REAL                                 :: mod_wind, buf
324
325!****************************************************************************************
326! Calculate the surface flux
327!
328!****************************************************************************************
329    DO i=1,knon
330       mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
331       buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
332       flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
333       flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
334    END DO
335
336  END SUBROUTINE calcul_flux_wind
337!
338!****************************************************************************************
339!
340END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.