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

Last change on this file since 5153 was 5153, checked in by abarral, 3 months ago

Revert FCTTRE to INCLUDE to assess impact of inlining
  • 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.2 KB
Line 
1! $Id: calcul_fluxs_mod.F90 5153 2024-07-31 16:20:03Z abarral $
2
3MODULE calcul_fluxs_mod
4
5  IMPLICIT NONE
6
7CONTAINS
8  SUBROUTINE calcul_fluxs(knon, nisurf, dtime, &
9          tsurf, p1lay, cal, beta, cdragh, cdragq, ps, &
10          precip_rain, precip_snow, snow, qsurf, &
11          radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, &
12          fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, &
13          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
14          sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
15
16    USE dimphy, ONLY: klon
17    USE indice_sol_mod
18    USE sens_heat_rain_m, ONLY: sens_heat_rain
19    USE lmdz_clesphys
20    USE lmdz_yoethf
21
22    USE lmdz_yomcst
23
24    IMPLICIT NONE
25 INCLUDE "FCTTRE.h"
26
27    ! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
28    ! une temperature de surface (au cas ou ok_veget = false)
29
30    ! L. Fairhead 4/2000
31
32    ! input:
33    !   knon         nombre de points a traiter
34    !   nisurf       surface a traiter
35    !   tsurf        temperature de surface
36    !   p1lay        pression 1er niveau (milieu de couche)
37    !   cal          capacite calorifique du sol
38    !   beta         evap reelle
39    !   cdragh       coefficient d'echange temperature
40    !   cdragq       coefficient d'echange evaporation
41    !   ps           pression au sol
42    !   precip_rain  precipitations liquides
43    !   precip_snow  precipitations solides
44    !   snow         champs hauteur de neige
45    !   runoff       runoff en cas de trop plein
46    !   petAcoef     coeff. A de la resolution de la CL pour t
47    !   peqAcoef     coeff. A de la resolution de la CL pour q
48    !   petBcoef     coeff. B de la resolution de la CL pour t
49    !   peqBcoef     coeff. B de la resolution de la CL pour q
50    !   radsol       rayonnement net aus sol (LW + SW)
51    !   dif_grnd     coeff. diffusion vers le sol profond
52
53    ! output:
54    !   tsurf_new    temperature au sol
55    !   qsurf        humidite de l'air au dessus du sol
56    !   fluxsens     flux de chaleur sensible
57    !   fluxlat      flux de chaleur latente
58    !   dflux_s      derivee du flux de chaleur sensible / Ts
59    !   dflux_l      derivee du flux de chaleur latente  / Ts
60    !   sens_prec_liq flux sensible lié aux echanges de precipitations liquides
61    !   sens_prec_sol                                    precipitations solides
62    !   lat_prec_liq  flux latent lié aux echanges de precipitations liquides
63    !   lat_prec_sol                                  precipitations solides
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
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    REAL, INTENT(IN), optional :: rhoa(:) ! (knon)
79    ! density of moist air  (kg / m3)
80
81    ! Parametres entree-sorties
82    !****************************************************************************************
83    REAL, DIMENSION(klon), INTENT(INOUT) :: snow  ! snow pas utile
84
85    ! Parametres sorties
86    !****************************************************************************************
87    REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
88    REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new, evap, fluxsens, fluxlat
89    REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
90    REAL, INTENT(OUT), OPTIONAL :: sens_prec_liq(:), sens_prec_sol(:) ! (knon)
91    REAL, DIMENSION(klon), OPTIONAL :: lat_prec_liq, lat_prec_sol
92
93    ! Variables locales
94    !****************************************************************************************
95    INTEGER :: i
96    REAL, DIMENSION(klon) :: zx_mh, zx_nh, zx_oh
97    REAL, DIMENSION(klon) :: zx_mq, zx_nq, zx_oq
98    REAL, DIMENSION(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat
99    REAL, DIMENSION(klon) :: zx_sl, zx_coefh, zx_coefq, zx_wind
100    REAL, DIMENSION(klon) :: d_ts
101    REAL :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
102    REAL :: qsat_new, q1_new
103    REAL, PARAMETER :: t_grnd = 271.35, t_coup = 273.15
104    REAL, PARAMETER :: max_eau_sol = 150.0
105    CHARACTER (len = 20) :: modname = 'calcul_fluxs'
106    LOGICAL :: fonte_neige
107    LOGICAL, SAVE :: check = .FALSE.
108    !$OMP THREADPRIVATE(check)
109
110    ! End definition
111    !****************************************************************************************
112
113    IF (check) WRITE(*, *)'Entree ', modname, ' surface = ', nisurf
114
115    IF (check) THEN
116      WRITE(*, *)' radsol (min, max)', &
117              MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
118    ENDIF
119
120    ! Traitement neige et humidite du sol
121    !****************************************************************************************
122
123    !!$  WRITE(*,*)'test calcul_flux, surface ', nisurf
124    !!PB test
125    !!$    if (nisurf == is_oce) THEN
126    !!$      snow = 0.
127    !!$      qsol = max_eau_sol
128    !!$    else
129    !!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
130    !!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
131    !!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)
132    !!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
133    !!$    endif
134    !!$    IF (nisurf /= is_ter) qsol = max_eau_sol
135
136    ! Initialisation
137    !****************************************************************************************
138    evap = 0.
139    fluxsens = 0.
140    fluxlat = 0.
141    dflux_s = 0.
142    dflux_l = 0.
143    IF (PRESENT(lat_prec_liq)) lat_prec_liq = 0.
144    IF (PRESENT(lat_prec_sol)) lat_prec_sol = 0.
145
146    ! zx_qs = qsat en kg/kg
147    !****************************************************************************************
148    DO i = 1, knon
149      zx_pkh(i) = (ps(i) / ps(i))**RKAPPA
150      IF (thermcep) THEN
151        zdelta = MAX(0., SIGN(1., rtt - tsurf(i)))
152        zcvm5 = R5LES * RLVTT * (1. - zdelta) + R5IES * RLSTT * zdelta
153        zcvm5 = zcvm5 / RCPD / (1.0 + RVTMP2 * q1lay(i))
154        zx_qs = r2es * FOEEW(tsurf(i), zdelta) / ps(i)
155        zx_qs = MIN(0.5, zx_qs)
156        zcor = 1. / (1. - retv * zx_qs)
157        zx_qs = zx_qs * zcor
158        zx_dq_s_dh = FOEDE(tsurf(i), zdelta, zcvm5, zx_qs, zcor) &
159                / RLVTT / zx_pkh(i)
160      ELSE
161        IF (tsurf(i)<t_coup) THEN
162          zx_qs = qsats(tsurf(i)) / ps(i)
163          zx_dq_s_dh = dqsats(tsurf(i), zx_qs) / RLVTT &
164                  / zx_pkh(i)
165        ELSE
166          zx_qs = qsatl(tsurf(i)) / ps(i)
167          zx_dq_s_dh = dqsatl(tsurf(i), zx_qs) / RLVTT &
168                  / zx_pkh(i)
169        ENDIF
170      ENDIF
171      zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
172      zx_qsat(i) = zx_qs
173      zx_wind(i) = min_wind_speed + SQRT(gustiness(i) + u1lay(i)**2 + v1lay(i)**2)
174      zx_coefh(i) = cdragh(i) * zx_wind(i) * p1lay(i) / (RD * t1lay(i))
175      zx_coefq(i) = cdragq(i) * zx_wind(i) * p1lay(i) / (RD * t1lay(i))
176      !      zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)**2+v1lay(i)**2) &
177      !                * p1lay(i)/(RD*t1lay(i))
178      !      zx_coefh(i) = cdragh(i) * zx_wind(i)
179      !      zx_coefq(i) = cdragq(i) * zx_wind(i)
180    ENDDO
181
182
183    ! === Calcul de la temperature de surface ===
184    ! zx_sl = chaleur latente d'evaporation ou de sublimation
185    !****************************************************************************************
186
187    DO i = 1, knon
188      zx_sl(i) = RLVTT
189      IF (tsurf(i) < RTT) zx_sl(i) = RLSTT
190    ENDDO
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      ! Y'a-t-il fonte de neige?
218
219      !    fonte_neige = (nisurf /= is_oce) .AND. &
220      !     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
221      !     & .AND. (tsurf_new(i) >= RTT)
222      !    if (fonte_neige) tsurf_new(i) = RTT
223      d_ts(i) = tsurf_new(i) - tsurf(i)
224      !    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
225      !    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
226
227      !== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
228      !== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
229      evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
230      fluxlat(i) = - evap(i) * zx_sl(i)
231      fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
232
233      ! Derives des flux dF/dTs (W m-2 K-1):
234      dflux_s(i) = zx_nh(i)
235      dflux_l(i) = (zx_sl(i) * zx_nq(i))
236
237      ! Nouvelle valeure de l'humidite au dessus du sol
238      qsat_new = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
239      q1_new = peqAcoef(i) - peqBcoef(i) * evap(i) * dtime
240      qsurf(i) = q1_new * (1. - beta(i)) + beta(i) * qsat_new
241
242      ! en cas de fonte de neige
243
244      !    if (fonte_neige) THEN
245      !      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
246      !     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
247      !     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
248      !      bilan_f = max(0., bilan_f)
249      !      fq_fonte = bilan_f / zx_sl(i)
250      !      snow(i) = max(0., snow(i) - fq_fonte * dtime)
251      !      qsol(i) = qsol(i) + (fq_fonte * dtime)
252      !    endif
253      !!$    if (nisurf == is_ter)  &
254      !!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
255      !!$    qsol(i) = min(qsol(i), max_eau_sol)
256
257      ! calcul de l'enthalpie des precipitations liquides et solides
258      IF (PRESENT(sens_prec_liq)) sens_prec_liq(i) &
259              = - sens_heat_rain(precip_rain(i) + precip_snow(i), t1lay(i), &
260                      q1lay(i), rhoa(i), rlvtt, tsurf_new(i), ps(i))
261      IF (PRESENT(sens_prec_sol)) sens_prec_sol(i) = 0.
262      ! On calcule par rapport a T=0
263      !! sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * precip_rain(i)
264      !! sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * precip_snow(i)
265
266      IF (PRESENT(lat_prec_liq))                    &
267              lat_prec_liq(i) = precip_rain(i) * (RLVTT - RLVTT)
268      IF (PRESENT(lat_prec_sol))                    &
269              lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
270    ENDDO
271
272
273    !**************************************************************************
274
275  END SUBROUTINE calcul_fluxs
276
277  !****************************************************************************************
278
279  SUBROUTINE calcul_flux_wind(knon, dtime, &
280          u0, v0, u1, v1, gustiness, cdrag_m, &
281          AcoefU, AcoefV, BcoefU, BcoefV, &
282          p1lay, t1lay, &
283          flux_u1, flux_v1)
284
285    USE dimphy
286    USE lmdz_clesphys
287    USE lmdz_yomcst
288
289    IMPLICIT NONE
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.