source: LMDZ6/branches/contrails/libf/phylmd/calcul_fluxs_mod.f90 @ 5456

Last change on this file since 5456 was 5285, checked in by abarral, 2 months ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules
  • 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 5285 2024-10-28 13:33:29Z aborella $
3!
4MODULE calcul_fluxs_mod
5
6  USE clesphys_mod_h
7    IMPLICIT NONE
8
9CONTAINS
10  SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
11       tsurf, p1lay, cal, beta, cdragh, cdragq, ps, &
12       precip_rain, precip_snow, snow, qsurf, &
13       radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, &
14       fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, &
15       tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
16       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
17
18
19    USE dimphy, ONLY : klon
20    USE indice_sol_mod
21    use sens_heat_rain_m, only: sens_heat_rain
22    USE yomcst_mod_h
23    USE yoethf_mod_h
24
25
26! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
27! une temperature de surface (au cas ou ok_veget = false)
28!
29! L. Fairhead 4/2000
30!
31! input:
32!   knon         nombre de points a traiter
33!   nisurf       surface a traiter
34!   tsurf        temperature de surface
35!   p1lay        pression 1er niveau (milieu de couche)
36!   cal          capacite calorifique du sol
37!   beta         evap reelle
38!   cdragh       coefficient d'echange temperature
39!   cdragq       coefficient d'echange evaporation
40!   ps           pression au sol
41!   precip_rain  precipitations liquides
42!   precip_snow  precipitations solides
43!   snow         champs hauteur de neige
44!   runoff       runoff en cas de trop plein
45!   petAcoef     coeff. A de la resolution de la CL pour t
46!   peqAcoef     coeff. A de la resolution de la CL pour q
47!   petBcoef     coeff. B de la resolution de la CL pour t
48!   peqBcoef     coeff. B de la resolution de la CL pour q
49!   radsol       rayonnement net aus sol (LW + SW)
50!   dif_grnd     coeff. diffusion vers le sol profond
51!
52! output:
53!   tsurf_new    temperature au sol
54!   qsurf        humidite de l'air au dessus du sol
55!   fluxsens     flux de chaleur sensible
56!   fluxlat      flux de chaleur latente
57!   dflux_s      derivee du flux de chaleur sensible / Ts
58!   dflux_l      derivee du flux de chaleur latente  / Ts
59!   sens_prec_liq flux sensible li� aux echanges de precipitations liquides
60!   sens_prec_sol                                    precipitations solides
61!   lat_prec_liq  flux latent li� aux echanges de precipitations liquides
62!   lat_prec_sol                                  precipitations solides
63
64    INCLUDE "FCTTRE.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 clesphys_mod_h
291    USE yomcst_mod_h
292USE dimphy
293
294
295! Input arguments
296!****************************************************************************************
297    INTEGER, INTENT(IN)                  :: knon
298    REAL, INTENT(IN)                     :: dtime
299    REAL, DIMENSION(klon), INTENT(IN)    :: u0, v0  ! u and v at niveau 0
300    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness  ! u and v at niveau 1
301    REAL, DIMENSION(klon), INTENT(IN)    :: cdrag_m ! cdrag pour momentum
302    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
303    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay   ! pression 1er niveau (milieu de couche)
304    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay   ! temperature
305! Output arguments
306!****************************************************************************************
307    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1
308    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v1
309
310! Local variables
311!****************************************************************************************
312    INTEGER                              :: i
313    REAL                                 :: mod_wind, buf
314
315!****************************************************************************************
316! Calculate the surface flux
317!
318!****************************************************************************************
319    DO i=1,knon
320       mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
321       buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
322       flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
323       flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
324    END DO
325
326  END SUBROUTINE calcul_flux_wind
327!
328!****************************************************************************************
329!
330END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.