source: LMDZ6/trunk/libf/phylmd/calcul_fluxs_mod.f90 @ 5496

Last change on this file since 5496 was 5486, checked in by evignon, 5 days ago

inclusion d'un diagnostique de la sublimation de la glace sur les landice
pour la conservation de l'eau
  • 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: 12.9 KB
Line 
1!
2! $Id: calcul_fluxs_mod.f90 5486 2025-01-17 14:38:37Z evignon $
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    ENDDO
180
181
182! === Calcul de la temperature de surface ===
183! zx_sl = chaleur latente d'evaporation ou de sublimation
184!****************************************************************************************
185
186    DO i = 1, knon
187       zx_sl(i) = RLVTT
188       IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
189    ENDDO
190   
191
192    DO i = 1, knon
193! Q
194       zx_oq(i) = 1. - (beta(i) * zx_coefq(i) * peqBcoef(i) * dtime)
195       zx_mq(i) = beta(i) * zx_coefq(i) * &
196            (peqAcoef(i) -             &
197! conv num avec precedente version
198            fqsat * zx_qsat(i) + fqsat * zx_dq_s_dt(i) * tsurf(i))  &
199!           fqsat * ( zx_qsat(i) - zx_dq_s_dt(i) * tsurf(i)) ) &
200            / zx_oq(i)
201       zx_nq(i) = beta(i) * zx_coefq(i) * (- fqsat * zx_dq_s_dt(i)) &
202            / zx_oq(i)
203       
204! H
205       zx_oh(i) = 1. - (zx_coefh(i) * petBcoef(i) * dtime)
206       zx_mh(i) = zx_coefh(i) * petAcoef(i) / zx_oh(i)
207       zx_nh(i) = - (zx_coefh(i) * RCPD * zx_pkh(i))/ zx_oh(i)
208     
209! Tsurface
210       tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
211            (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
212            + dif_grnd(i) * t_grnd * dtime)/ &
213            ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
214            zx_nh(i) + zx_sl(i) * zx_nq(i)) & 
215            + dtime * dif_grnd(i))
216
217!
218! Y'a-t-il fonte de neige?
219!
220!    fonte_neige = (nisurf /= is_oce) .AND. &
221!     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
222!     & .AND. (tsurf_new(i) >= RTT)
223!    if (fonte_neige) tsurf_new(i) = RTT 
224       d_ts(i) = tsurf_new(i) - tsurf(i)
225!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
226!    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
227
228!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
229!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
230       evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
231       fluxlat(i) = - evap(i) * zx_sl(i)
232       fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
233       
234! Derives des flux dF/dTs (W m-2 K-1):
235       dflux_s(i) = zx_nh(i)
236       dflux_l(i) = (zx_sl(i) * zx_nq(i))
237
238! Nouvelle valeure de l'humidite au dessus du sol
239       qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
240       q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
241       qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
242!
243! en cas de fonte de neige
244!
245!    if (fonte_neige) then
246!      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
247!     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
248!     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
249!      bilan_f = max(0., bilan_f)
250!      fq_fonte = bilan_f / zx_sl(i)
251!      snow(i) = max(0., snow(i) - fq_fonte * dtime)
252!      qsol(i) = qsol(i) + (fq_fonte * dtime)
253!    endif
254!!$    if (nisurf == is_ter)  &
255!!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
256!!$    qsol(i) = min(qsol(i), max_eau_sol)
257!
258! calcul de l'enthalpie des precipitations liquides et solides
259       if (PRESENT(sens_prec_liq)) sens_prec_liq(i) &
260            = - sens_heat_rain(precip_rain(i) + precip_snow(i), t1lay(i), &
261            q1lay(i), rhoa(i), rlvtt, tsurf_new(i), ps(i))
262       if (PRESENT(sens_prec_sol)) sens_prec_sol(i) = 0.
263       ! On calcule par rapport a T=0
264       !! sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * precip_rain(i)
265       !! sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * precip_snow(i)
266
267       if (PRESENT(lat_prec_liq))                    &
268         lat_prec_liq(i) =  precip_rain(i) * (RLVTT - RLVTT)
269       if (PRESENT(lat_prec_sol))                    &
270         lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
271    ENDDO
272
273   
274!**************************************************************************
275!
276  END SUBROUTINE calcul_fluxs
277!
278!****************************************************************************************
279!
280  SUBROUTINE calcul_flux_wind(knon, dtime, &
281       u0, v0, u1, v1, gustiness, cdrag_m, &
282       AcoefU, AcoefV, BcoefU, BcoefV, &
283       p1lay, t1lay, &
284       flux_u1, flux_v1)
285
286    USE clesphys_mod_h
287    USE yomcst_mod_h
288USE dimphy
289
290
291! Input arguments
292!****************************************************************************************
293    INTEGER, INTENT(IN)                  :: knon
294    REAL, INTENT(IN)                     :: dtime
295    REAL, DIMENSION(klon), INTENT(IN)    :: u0, v0  ! u and v at niveau 0
296    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness  ! u and v at niveau 1
297    REAL, DIMENSION(klon), INTENT(IN)    :: cdrag_m ! cdrag pour momentum
298    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
299    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay   ! pression 1er niveau (milieu de couche)
300    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay   ! temperature
301! Output arguments
302!****************************************************************************************
303    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1
304    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v1
305
306! Local variables
307!****************************************************************************************
308    INTEGER                              :: i
309    REAL                                 :: mod_wind, buf
310
311!****************************************************************************************
312! Calculate the surface flux
313!
314!****************************************************************************************
315    DO i=1,knon
316       mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
317       buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
318       flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
319       flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
320    END DO
321
322  END SUBROUTINE calcul_flux_wind
323!
324!****************************************************************************************
325!
326END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.