source: LMDZ6/trunk/libf/phylmd/calcul_fluxs_mod.F90 @ 3409

Last change on this file since 3409 was 3319, checked in by jyg, 7 years ago

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