source: lmdz_wrf/trunk/WRFV3/lmdz/calcul_fluxs_mod.F90 @ 2294

Last change on this file since 2294 was 183, checked in by lfita, 10 years ago

Removing check pritings
Removing nullification of MixingRatioValues?(:,:,2) = 0 and MixingTendRatioValues?(:,:,2)
NO plul related to te absence of oliq at the initial conditions. If one introduce fake ones it works!
  • Property svn:executable set to *
File size: 11.0 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    USE indice_sol_mod
15
16! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
17! une temperature de surface (au cas ou ok_veget = false)
18!
19! L. Fairhead 4/2000
20!
21! input:
22!   knon         nombre de points a traiter
23!   nisurf       surface a traiter
24!   tsurf        temperature de surface
25!   p1lay        pression 1er niveau (milieu de couche)
26!   cal          capacite calorifique du sol
27!   beta         evap reelle
28!   coef1lay     coefficient d'echange
29!   ps           pression au sol
30!   precip_rain  precipitations liquides
31!   precip_snow  precipitations solides
32!   snow         champs hauteur de neige
33!   runoff       runoff en cas de trop plein
34!   petAcoef     coeff. A de la resolution de la CL pour t
35!   peqAcoef     coeff. A de la resolution de la CL pour q
36!   petBcoef     coeff. B de la resolution de la CL pour t
37!   peqBcoef     coeff. B de la resolution de la CL pour q
38!   radsol       rayonnement net aus sol (LW + SW)
39!   dif_grnd     coeff. diffusion vers le sol profond
40!
41! output:
42!   tsurf_new    temperature au sol
43!   qsurf        humidite de l'air au dessus du sol
44!   fluxsens     flux de chaleur sensible
45!   fluxlat      flux de chaleur latente
46!   dflux_s      derivee du flux de chaleur sensible / Ts
47!   dflux_l      derivee du flux de chaleur latente  / Ts
48!
49
50    INCLUDE "YOETHF.h"
51    INCLUDE "FCTTRE.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!Lluis
96    INTEGER                                              :: lpt
97    lpt = 550
98
99    IF (check) WRITE(*,*)'Entree ', modname,' surface = ',nisurf
100   
101    IF (check) THEN
102       WRITE(*,*)' radsol (min, max)', &
103            MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
104    ENDIF
105 
106! Traitement neige et humidite du sol
107!****************************************************************************************
108!
109!!$  WRITE(*,*)'test calcul_flux, surface ', nisurf
110!!PB test
111!!$    if (nisurf == is_oce) then
112!!$      snow = 0.
113!!$      qsol = max_eau_sol
114!!$    else
115!!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
116!!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
117!!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)
118!!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
119!!$    endif
120!!$    IF (nisurf /= is_ter) qsol = max_eau_sol
121
122
123!
124! Initialisation
125!****************************************************************************************
126    evap = 0.
127    fluxsens=0.
128    fluxlat=0.
129    dflux_s = 0.
130    dflux_l = 0.       
131!
132! zx_qs = qsat en kg/kg
133!****************************************************************************************
134    DO i = 1, knon
135       zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
136       IF (thermcep) THEN
137          zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
138          zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
139          zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
140          zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
141          zx_qs=MIN(0.5,zx_qs)
142          zcor=1./(1.-retv*zx_qs)
143          zx_qs=zx_qs*zcor
144          zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
145               /RLVTT / zx_pkh(i)
146       ELSE
147          IF (tsurf(i).LT.t_coup) THEN
148             zx_qs = qsats(tsurf(i)) / ps(i)
149             zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
150                  / zx_pkh(i)
151          ELSE
152             zx_qs = qsatl(tsurf(i)) / ps(i)
153             zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
154                  / zx_pkh(i)
155          ENDIF
156       ENDIF
157       zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
158       zx_qsat(i) = zx_qs
159       zx_coef(i) = coef1lay(i) * &
160            (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) * &
161            p1lay(i)/(RD*t1lay(i))
162       
163    ENDDO
164
165
166! === Calcul de la temperature de surface ===
167! zx_sl = chaleur latente d'evaporation ou de sublimation
168!****************************************************************************************
169
170    DO i = 1, knon
171       zx_sl(i) = RLVTT
172       IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
173       zx_k1(i) = zx_coef(i)
174    ENDDO
175   
176
177    DO i = 1, knon
178! Q
179       zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
180       zx_mq(i) = beta(i) * zx_k1(i) * &
181            (peqAcoef(i) - zx_qsat(i) + &
182            zx_dq_s_dt(i) * tsurf(i)) &
183            / zx_oq(i)
184       zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
185            / zx_oq(i)
186       
187! H
188       zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
189       zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
190       zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
191     
192! Tsurface
193       tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
194            (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
195            + dif_grnd(i) * t_grnd * dtime)/ &
196            ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
197            zx_nh(i) + zx_sl(i) * zx_nq(i)) & 
198            + dtime * dif_grnd(i))
199
200!
201! Y'a-t-il fonte de neige?
202!
203!    fonte_neige = (nisurf /= is_oce) .AND. &
204!     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
205!     & .AND. (tsurf_new(i) >= RTT)
206!    if (fonte_neige) tsurf_new(i) = RTT 
207       d_ts(i) = tsurf_new(i) - tsurf(i)
208!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
209!    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
210
211!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
212!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
213       evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
214       fluxlat(i) = - evap(i) * zx_sl(i)
215       fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
216       
217! Derives des flux dF/dTs (W m-2 K-1):
218       dflux_s(i) = zx_nh(i)
219       dflux_l(i) = (zx_sl(i) * zx_nq(i))
220
221! Nouvelle valeure de l'humidite au dessus du sol
222       qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
223       q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
224       qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
225!
226! en cas de fonte de neige
227!
228!    if (fonte_neige) then
229!      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
230!     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
231!     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
232!      bilan_f = max(0., bilan_f)
233!      fq_fonte = bilan_f / zx_sl(i)
234!      snow(i) = max(0., snow(i) - fq_fonte * dtime)
235!      qsol(i) = qsol(i) + (fq_fonte * dtime)
236!    endif
237!!$    if (nisurf == is_ter)  &
238!!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
239!!$    qsol(i) = min(qsol(i), max_eau_sol)
240    ENDDO
241!
242!****************************************************************************************
243!
244  END SUBROUTINE calcul_fluxs
245!
246!****************************************************************************************
247!
248  SUBROUTINE calcul_flux_wind(knon, dtime, &
249       u0, v0, u1, v1, cdrag_m, &
250       AcoefU, AcoefV, BcoefU, BcoefV, &
251       p1lay, t1lay, &
252       flux_u1, flux_v1)
253
254    USE dimphy
255    INCLUDE "YOMCST.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  ! 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 = 1.0 + SQRT((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.