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

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

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