source: LMDZ6/branches/Amaury_dev/libf/phylmd/calcul_fluxs_mod.F90 @ 5134

Last change on this file since 5134 was 5117, checked in by abarral, 4 months ago

rename modules properly lmdz_*
move some unused files to obsolete/
(lint) uppercase fortran keywords
  • 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 5117 2024-07-24 14:23:34Z abarral $
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! Initialisation
138!****************************************************************************************
139    evap = 0.
140    fluxsens=0.
141    fluxlat=0.
142    dflux_s = 0.
143    dflux_l = 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)<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) < 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! Y'a-t-il fonte de neige?
220
221!    fonte_neige = (nisurf /= is_oce) .AND. &
222!     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
223!     & .AND. (tsurf_new(i) >= RTT)
224!    if (fonte_neige) tsurf_new(i) = RTT 
225       d_ts(i) = tsurf_new(i) - tsurf(i)
226!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
227!    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
228
229!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
230!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
231       evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
232       fluxlat(i) = - evap(i) * zx_sl(i)
233       fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
234       
235! Derives des flux dF/dTs (W m-2 K-1):
236       dflux_s(i) = zx_nh(i)
237       dflux_l(i) = (zx_sl(i) * zx_nq(i))
238
239! Nouvelle valeure de l'humidite au dessus du sol
240       qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
241       q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
242       qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
243
244! en cas de fonte de neige
245
246!    if (fonte_neige) THEN
247!      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
248!     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
249!     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
250!      bilan_f = max(0., bilan_f)
251!      fq_fonte = bilan_f / zx_sl(i)
252!      snow(i) = max(0., snow(i) - fq_fonte * dtime)
253!      qsol(i) = qsol(i) + (fq_fonte * dtime)
254!    endif
255!!$    if (nisurf == is_ter)  &
256!!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
257!!$    qsol(i) = min(qsol(i), max_eau_sol)
258
259! calcul de l'enthalpie des precipitations liquides et solides
260       IF (PRESENT(sens_prec_liq)) sens_prec_liq(i) &
261            = - sens_heat_rain(precip_rain(i) + precip_snow(i), t1lay(i), &
262            q1lay(i), rhoa(i), rlvtt, tsurf_new(i), ps(i))
263       IF (PRESENT(sens_prec_sol)) sens_prec_sol(i) = 0.
264       ! On calcule par rapport a T=0
265       !! sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * precip_rain(i)
266       !! sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * precip_snow(i)
267
268       IF (PRESENT(lat_prec_liq))                    &
269         lat_prec_liq(i) =  precip_rain(i) * (RLVTT - RLVTT)
270       IF (PRESENT(lat_prec_sol))                    &
271         lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
272    ENDDO
273
274   
275!**************************************************************************
276
277  END SUBROUTINE calcul_fluxs
278
279!****************************************************************************************
280
281  SUBROUTINE calcul_flux_wind(knon, dtime, &
282       u0, v0, u1, v1, gustiness, cdrag_m, &
283       AcoefU, AcoefV, BcoefU, BcoefV, &
284       p1lay, t1lay, &
285       flux_u1, flux_v1)
286
287    USE dimphy
288    INCLUDE "YOMCST.h"
289    INCLUDE "clesphys.h"
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.