source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/calcul_fluxs_mod.F90 @ 5444

Last change on this file since 5444 was 1106, checked in by lguez, 16 years ago

"comconst.h", and "comgeom2.h" are now both fixed and free form.
Removed calls to procedure "flush".
Corrected kinds of constants which appeared as arguments to "min" or
"max" (all arguments are now of the same type and kind).
  • 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    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!
289END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.