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

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

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