source: LMDZ5/trunk/libf/phylmd/calcul_fluxs_mod.F90 @ 2245

Last change on this file since 2245 was 2240, checked in by fhourdin, 9 years ago

Revisite de la formule des flux de surface
(en priorité sur l'océan) en tenant compte des bourrasques de
vent et de la différence entre les hauteurs de rugosités pour
la quantité de mouvement, l'enthalpie et éventuellement l'humidité.

Etape 1 :
Introduction d'un calcul de gustiness dans la physique
gustiness(:)=f_gust_bl * ale_bl + f_gust_wk * ame_wk
Cette variable est passée ensuite jusqu'au fin fond de la couche limite.
L'étape 1 est prête à commettre, ne nécessite pas de nouvelles
variables dans les startphy et assure la convergence numérique.

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