source: LMDZ4/trunk/libf/phylmd/calcul_fluxs_mod.F90 @ 1081

Last change on this file since 1081 was 1067, checked in by Laurent Fairhead, 16 years ago
  • Modifications lie au premier niveau du modele pour la diffusion turbulent

du vent.

  • Preparation pour un couplage des courrant oceaniques.

JG
  • 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
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, &
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       CALL flush(6)
102    ENDIF
103 
104! Traitement neige et humidite du sol
105!****************************************************************************************
106!
107!!$  WRITE(*,*)'test calcul_flux, surface ', nisurf
108!!PB test
109!!$    if (nisurf == is_oce) then
110!!$      snow = 0.
111!!$      qsol = max_eau_sol
112!!$    else
113!!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
114!!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
115!!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)
116!!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
117!!$    endif
118!!$    IF (nisurf /= is_ter) qsol = max_eau_sol
119
120
121!
122! Initialisation
123!****************************************************************************************
124    evap = 0.
125    fluxsens=0.
126    fluxlat=0.
127    dflux_s = 0.
128    dflux_l = 0.       
129!
130! zx_qs = qsat en kg/kg
131!****************************************************************************************
132    DO i = 1, knon
133       zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
134       IF (thermcep) THEN
135          zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
136          zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
137          zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
138          zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
139          zx_qs=MIN(0.5,zx_qs)
140          zcor=1./(1.-retv*zx_qs)
141          zx_qs=zx_qs*zcor
142          zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
143               /RLVTT / zx_pkh(i)
144       ELSE
145          IF (tsurf(i).LT.t_coup) THEN
146             zx_qs = qsats(tsurf(i)) / ps(i)
147             zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
148                  / zx_pkh(i)
149          ELSE
150             zx_qs = qsatl(tsurf(i)) / ps(i)
151             zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
152                  / zx_pkh(i)
153          ENDIF
154       ENDIF
155       zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
156       zx_qsat(i) = zx_qs
157       zx_coef(i) = coef1lay(i) * &
158            (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) * &
159            p1lay(i)/(RD*t1lay(i))
160       
161    ENDDO
162
163
164! === Calcul de la temperature de surface ===
165! zx_sl = chaleur latente d'evaporation ou de sublimation
166!****************************************************************************************
167
168    DO i = 1, knon
169       zx_sl(i) = RLVTT
170       IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
171       zx_k1(i) = zx_coef(i)
172    ENDDO
173   
174
175    DO i = 1, knon
176! Q
177       zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
178       zx_mq(i) = beta(i) * zx_k1(i) * &
179            (peqAcoef(i) - zx_qsat(i) + &
180            zx_dq_s_dt(i) * tsurf(i)) &
181            / zx_oq(i)
182       zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
183            / zx_oq(i)
184       
185! H
186       zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
187       zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
188       zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
189     
190! Tsurface
191       tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
192            (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
193            + dif_grnd(i) * t_grnd * dtime)/ &
194            ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
195            zx_nh(i) + zx_sl(i) * zx_nq(i)) & 
196            + dtime * dif_grnd(i))
197
198!
199! Y'a-t-il fonte de neige?
200!
201!    fonte_neige = (nisurf /= is_oce) .AND. &
202!     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
203!     & .AND. (tsurf_new(i) >= RTT)
204!    if (fonte_neige) tsurf_new(i) = RTT 
205       d_ts(i) = tsurf_new(i) - tsurf(i)
206!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
207!    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
208
209!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
210!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
211       evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
212       fluxlat(i) = - evap(i) * zx_sl(i)
213       fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
214       
215! Derives des flux dF/dTs (W m-2 K-1):
216       dflux_s(i) = zx_nh(i)
217       dflux_l(i) = (zx_sl(i) * zx_nq(i))
218
219! Nouvelle valeure de l'humidite au dessus du sol
220       qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
221       q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
222       qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
223!
224! en cas de fonte de neige
225!
226!    if (fonte_neige) then
227!      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
228!     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
229!     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
230!      bilan_f = max(0., bilan_f)
231!      fq_fonte = bilan_f / zx_sl(i)
232!      snow(i) = max(0., snow(i) - fq_fonte * dtime)
233!      qsol(i) = qsol(i) + (fq_fonte * dtime)
234!    endif
235!!$    if (nisurf == is_ter)  &
236!!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
237!!$    qsol(i) = min(qsol(i), max_eau_sol)
238    ENDDO
239!
240!****************************************************************************************
241!
242  END SUBROUTINE calcul_fluxs
243!
244!****************************************************************************************
245!
246  SUBROUTINE calcul_flux_wind(knon, dtime, &
247       u0, v0, u1, v1, cdrag_m, &
248       AcoefU, AcoefV, BcoefU, BcoefV, &
249       p1lay, t1lay, &
250       flux_u1, flux_v1)
251
252    USE dimphy
253    INCLUDE "YOMCST.h"
254
255! Input arguments
256!****************************************************************************************
257    INTEGER, INTENT(IN)                  :: knon
258    REAL, INTENT(IN)                     :: dtime
259    REAL, DIMENSION(klon), INTENT(IN)    :: u0, v0  ! u and v at niveau 0
260    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1  ! u and v at niveau 1
261    REAL, DIMENSION(klon), INTENT(IN)    :: cdrag_m ! cdrag pour momentum
262    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
263    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay   ! pression 1er niveau (milieu de couche)
264    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay   ! temperature
265! Output arguments
266!****************************************************************************************
267    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1
268    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v1
269
270! Local variables
271!****************************************************************************************
272    INTEGER                              :: i
273    REAL                                 :: mod_wind, buf
274
275!****************************************************************************************
276! Calculate the surface flux
277!
278!****************************************************************************************
279    DO i=1,knon
280       mod_wind = 1.0 + SQRT((u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
281       buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
282       flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
283       flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
284    END DO
285
286  END SUBROUTINE calcul_flux_wind
287!
288!****************************************************************************************
289!
290END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.