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

Last change on this file since 3817 was 3815, checked in by lguez, 4 years ago

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