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

Last change on this file since 2488 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
RevLine 
[782]1!
2MODULE calcul_fluxs_mod
3
4
5CONTAINS
6  SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
[2298]7       tsurf, p1lay, cal, beta, cdragh, cdragq, ps, &
[782]8       precip_rain, precip_snow, snow, qsurf, &
[2298]9       radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, &
10       fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, &
[782]11       tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
12   
13    USE dimphy, ONLY : klon
[1795]14    USE indice_sol_mod
[782]15
[2298]16    INCLUDE "clesphys.h"
17
[782]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
[2298]30!   cdragh       coefficient d'echange temperature
31!   cdragq       coefficient d'echange evaporation
[782]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
[793]53    INCLUDE "YOETHF.h"
54    INCLUDE "FCTTRE.h"
55    INCLUDE "YOMCST.h"
[782]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
[2298]64    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf, p1lay, cal, beta, cdragh,cdragq
[782]65    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow ! pas utiles
66    REAL, DIMENSION(klon), INTENT(IN)    :: radsol, dif_grnd
[2298]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)
[782]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
[2298]85    REAL, DIMENSION(klon)                :: zx_pkh, zx_dq_s_dt, zx_qsat
86    REAL, DIMENSION(klon)                :: zx_sl, zx_coefh, zx_coefq, zx_wind
[782]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.
[2298]131    dflux_l = 0.
[782]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
[2298]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)
[782]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
[2298]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)) ) &
[782]188            / zx_oq(i)
[2298]189       zx_nq(i) = beta(i) * zx_coefq(i) * (- fqsat * zx_dq_s_dt(i)) &
[782]190            / zx_oq(i)
191       
192! H
[2298]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)
[782]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
[1067]250!
251!****************************************************************************************
252!
253  SUBROUTINE calcul_flux_wind(knon, dtime, &
[2298]254       u0, v0, u1, v1, gustiness, cdrag_m, &
[1067]255       AcoefU, AcoefV, BcoefU, BcoefV, &
256       p1lay, t1lay, &
257       flux_u1, flux_v1)
[782]258
[1067]259    USE dimphy
260    INCLUDE "YOMCST.h"
[2298]261    INCLUDE "clesphys.h"
[1067]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
[2298]268    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness  ! u and v at niveau 1
[1067]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
[2298]288       mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
[1067]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!
[782]298END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.