source: LMDZ5/trunk/libf/phylmd/calcul_fluxs_mod.F90 @ 4867

Last change on this file since 4867 was 2538, checked in by Laurent Fairhead, 9 years ago

Computation of heat fluxes associated with solid and liquid precipitations
over ocean and seaice. Quantities are sent to the coupler
LF
  • 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:executable set to *
  • Property svn:keywords set to Id
File size: 13.7 KB
RevLine 
[782]1!
[2538]2! $Id: calcul_fluxs_mod.F90 2538 2016-06-03 14:12:16Z idelkadi $
3!
[782]4MODULE calcul_fluxs_mod
5
6
7CONTAINS
8  SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
[2254]9       tsurf, p1lay, cal, beta, cdragh, cdragq, ps, &
[782]10       precip_rain, precip_snow, snow, qsurf, &
[2240]11       radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, &
[2254]12       fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, &
[2538]13       tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
14       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
15 
[782]16   
17    USE dimphy, ONLY : klon
[1785]18    USE indice_sol_mod
[782]19
[2240]20    INCLUDE "clesphys.h"
21
[782]22! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
23! une temperature de surface (au cas ou ok_veget = false)
24!
25! L. Fairhead 4/2000
26!
27! input:
28!   knon         nombre de points a traiter
29!   nisurf       surface a traiter
30!   tsurf        temperature de surface
31!   p1lay        pression 1er niveau (milieu de couche)
32!   cal          capacite calorifique du sol
33!   beta         evap reelle
[2254]34!   cdragh       coefficient d'echange temperature
35!   cdragq       coefficient d'echange evaporation
[782]36!   ps           pression au sol
37!   precip_rain  precipitations liquides
38!   precip_snow  precipitations solides
39!   snow         champs hauteur de neige
40!   runoff       runoff en cas de trop plein
41!   petAcoef     coeff. A de la resolution de la CL pour t
42!   peqAcoef     coeff. A de la resolution de la CL pour q
43!   petBcoef     coeff. B de la resolution de la CL pour t
44!   peqBcoef     coeff. B de la resolution de la CL pour q
45!   radsol       rayonnement net aus sol (LW + SW)
46!   dif_grnd     coeff. diffusion vers le sol profond
47!
48! output:
49!   tsurf_new    temperature au sol
50!   qsurf        humidite de l'air au dessus du sol
51!   fluxsens     flux de chaleur sensible
52!   fluxlat      flux de chaleur latente
53!   dflux_s      derivee du flux de chaleur sensible / Ts
54!   dflux_l      derivee du flux de chaleur latente  / Ts
[2538]55!   sens_prec_liq flux sensible lié aux echanges de precipitations liquides
56!   sens_prec_sol                                    precipitations solides
57!   lat_prec_liq  flux latent lié aux echanges de precipitations liquides
58!   lat_prec_sol                                  precipitations solides
[782]59
[793]60    INCLUDE "YOETHF.h"
61    INCLUDE "FCTTRE.h"
62    INCLUDE "YOMCST.h"
[782]63
64! Parametres d'entree
65!****************************************************************************************
66    INTEGER, INTENT(IN)                  :: knon, nisurf
67    REAL   , INTENT(IN)                  :: dtime
68    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
69    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
70    REAL, DIMENSION(klon), INTENT(IN)    :: ps, q1lay
[2254]71    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf, p1lay, cal, beta, cdragh,cdragq
[782]72    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow ! pas utiles
73    REAL, DIMENSION(klon), INTENT(IN)    :: radsol, dif_grnd
[2240]74    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay, u1lay, v1lay,gustiness
[2254]75    REAL,                  INTENT(IN)    :: fqsat ! correction factor on qsat (generally 0.98 over salty water, 1 everywhere else)
[782]76
77! Parametres entree-sorties
78!****************************************************************************************
79    REAL, DIMENSION(klon), INTENT(INOUT) :: snow  ! snow pas utile
80
81! Parametres sorties
82!****************************************************************************************
83    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
84    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new, evap, fluxsens, fluxlat
85    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l
[2538]86    REAL, DIMENSION(klon), OPTIONAL      :: sens_prec_liq, sens_prec_sol
87    REAL, DIMENSION(klon), OPTIONAL      :: lat_prec_liq, lat_prec_sol
[782]88
89! Variables locales
90!****************************************************************************************
91    INTEGER                              :: i
92    REAL, DIMENSION(klon)                :: zx_mh, zx_nh, zx_oh
93    REAL, DIMENSION(klon)                :: zx_mq, zx_nq, zx_oq
[2254]94    REAL, DIMENSION(klon)                :: zx_pkh, zx_dq_s_dt, zx_qsat
95    REAL, DIMENSION(klon)                :: zx_sl, zx_coefh, zx_coefq, zx_wind
[782]96    REAL, DIMENSION(klon)                :: d_ts
97    REAL                                 :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
98    REAL                                 :: qsat_new, q1_new
99    REAL, PARAMETER                      :: t_grnd = 271.35, t_coup = 273.15
100    REAL, PARAMETER                      :: max_eau_sol = 150.0
101    CHARACTER (len = 20)                 :: modname = 'calcul_fluxs'
102    LOGICAL                              :: fonte_neige
103    LOGICAL, SAVE                        :: check = .FALSE.
104    !$OMP THREADPRIVATE(check)
105
106! End definition
107!****************************************************************************************
108
109    IF (check) WRITE(*,*)'Entree ', modname,' surface = ',nisurf
110   
111    IF (check) THEN
112       WRITE(*,*)' radsol (min, max)', &
113            MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
114    ENDIF
115 
116! Traitement neige et humidite du sol
117!****************************************************************************************
118!
119!!$  WRITE(*,*)'test calcul_flux, surface ', nisurf
120!!PB test
121!!$    if (nisurf == is_oce) then
122!!$      snow = 0.
123!!$      qsol = max_eau_sol
124!!$    else
125!!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
126!!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
127!!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)
128!!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
129!!$    endif
130!!$    IF (nisurf /= is_ter) qsol = max_eau_sol
131
132
133!
134! Initialisation
135!****************************************************************************************
136    evap = 0.
137    fluxsens=0.
138    fluxlat=0.
139    dflux_s = 0.
[2254]140    dflux_l = 0.
[2538]141    if (PRESENT(sens_prec_liq)) sens_prec_liq = 0.
142    if (PRESENT(sens_prec_sol)) sens_prec_sol = 0.
143    if (PRESENT(lat_prec_liq)) lat_prec_liq = 0.
144    if (PRESENT(lat_prec_sol)) lat_prec_sol = 0.
[782]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).LT.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
[2254]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)
[782]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) .LT. RTT) zx_sl(i) = RLSTT
190    ENDDO
191   
192
193    DO i = 1, knon
194! Q
[2254]195       zx_oq(i) = 1. - (beta(i) * zx_coefq(i) * peqBcoef(i) * dtime)
196       zx_mq(i) = beta(i) * zx_coefq(i) * &
197            (peqAcoef(i) -             &
198! conv num avec precedente version
199            fqsat * zx_qsat(i) + fqsat * zx_dq_s_dt(i) * tsurf(i))  &
200!           fqsat * ( zx_qsat(i) - zx_dq_s_dt(i) * tsurf(i)) ) &
[782]201            / zx_oq(i)
[2254]202       zx_nq(i) = beta(i) * zx_coefq(i) * (- fqsat * zx_dq_s_dt(i)) &
[782]203            / zx_oq(i)
204       
205! H
[2254]206       zx_oh(i) = 1. - (zx_coefh(i) * petBcoef(i) * dtime)
207       zx_mh(i) = zx_coefh(i) * petAcoef(i) / zx_oh(i)
208       zx_nh(i) = - (zx_coefh(i) * RCPD * zx_pkh(i))/ zx_oh(i)
[782]209     
210! Tsurface
211       tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
212            (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
213            + dif_grnd(i) * t_grnd * dtime)/ &
214            ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
215            zx_nh(i) + zx_sl(i) * zx_nq(i)) & 
216            + dtime * dif_grnd(i))
217
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)
[2538]258!
259! calcul de l'enthalpie des precipitations liquides et solides
260!
261!       if (PRESENT(enth_prec_liq))                   &
262!       enth_prec_liq(i) = rcw * (t1lay(i) - tsurf(i)) * &
263!                          precip_rain(i)
264!       if (PRESENT(enth_prec_sol))                  &
265!       enth_prec_sol(i) = rcs * (t1lay(i) - tsurf(i)) * &
266!                          precip_snow(i)
267! On calcule par rapport a T=0
268       if (PRESENT(sens_prec_liq))                   &
269         sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * &
270                          precip_rain(i)
271       if (PRESENT(sens_prec_sol))                   &
272         sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * &
273                          precip_snow(i)
274       if (PRESENT(lat_prec_liq))                    &
275         lat_prec_liq(i) =  precip_rain(i) * (RLVTT - RLVTT)
276       if (PRESENT(lat_prec_sol))                    &
277         lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
[782]278    ENDDO
[2538]279
280   
281!       if (PRESENT(sens_prec_liq))                  &
282!          WRITE(*,*)' calculs_fluxs sens_prec_liq (min, max)', &
283!          MINVAL(sens_prec_liq(1:knon)), MAXVAL(sens_prec_liq(1:knon))
284!       if (PRESENT(sens_prec_sol))                  &
285!          WRITE(*,*)' calculs_fluxs sens_prec_sol (min, max)', &
286!          MINVAL(sens_prec_sol(1:knon)), MAXVAL(sens_prec_sol(1:knon))
287 
[782]288!
289!****************************************************************************************
290!
291  END SUBROUTINE calcul_fluxs
[1067]292!
293!****************************************************************************************
294!
295  SUBROUTINE calcul_flux_wind(knon, dtime, &
[2240]296       u0, v0, u1, v1, gustiness, cdrag_m, &
[1067]297       AcoefU, AcoefV, BcoefU, BcoefV, &
298       p1lay, t1lay, &
299       flux_u1, flux_v1)
[782]300
[1067]301    USE dimphy
302    INCLUDE "YOMCST.h"
[2240]303    INCLUDE "clesphys.h"
[1067]304
305! Input arguments
306!****************************************************************************************
307    INTEGER, INTENT(IN)                  :: knon
308    REAL, INTENT(IN)                     :: dtime
309    REAL, DIMENSION(klon), INTENT(IN)    :: u0, v0  ! u and v at niveau 0
[2240]310    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness  ! u and v at niveau 1
[1067]311    REAL, DIMENSION(klon), INTENT(IN)    :: cdrag_m ! cdrag pour momentum
312    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
313    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay   ! pression 1er niveau (milieu de couche)
314    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay   ! temperature
315! Output arguments
316!****************************************************************************************
317    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1
318    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v1
319
320! Local variables
321!****************************************************************************************
322    INTEGER                              :: i
323    REAL                                 :: mod_wind, buf
324
325!****************************************************************************************
326! Calculate the surface flux
327!
328!****************************************************************************************
329    DO i=1,knon
[2240]330       mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
[1067]331       buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
332       flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
333       flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
334    END DO
335
336  END SUBROUTINE calcul_flux_wind
337!
338!****************************************************************************************
339!
[782]340END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.