source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phylmd/calcul_fluxs_mod.F90 @ 5360

Last change on this file since 5360 was 3331, checked in by acozic, 7 years ago

Add modification for isotopes

  • 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
Line 
1!
2! $Id: calcul_fluxs_mod.F90 3331 2018-05-22 14:30:21Z abarral $
3!
4MODULE calcul_fluxs_mod
5
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)
15 
16   
17    USE dimphy, ONLY : klon
18    USE indice_sol_mod
19
20    INCLUDE "clesphys.h"
21
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
34!   cdragh       coefficient d'echange temperature
35!   cdragq       coefficient d'echange evaporation
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
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
59
60    INCLUDE "YOETHF.h"
61    INCLUDE "FCTTRE.h"
62    INCLUDE "YOMCST.h"
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
71    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf, p1lay, cal, beta, cdragh,cdragq
72    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow ! pas utiles
73    REAL, DIMENSION(klon), INTENT(IN)    :: radsol, dif_grnd
74    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay, u1lay, v1lay,gustiness
75    REAL,                  INTENT(IN)    :: fqsat ! correction factor on qsat (generally 0.98 over salty water, 1 everywhere else)
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
86    REAL, DIMENSION(klon), OPTIONAL      :: sens_prec_liq, sens_prec_sol
87    REAL, DIMENSION(klon), OPTIONAL      :: lat_prec_liq, lat_prec_sol
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
94    REAL, DIMENSION(klon)                :: zx_pkh, zx_dq_s_dt, zx_qsat
95    REAL, DIMENSION(klon)                :: zx_sl, zx_coefh, zx_coefq, zx_wind
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.
140    dflux_l = 0.
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.
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
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) .LT. RTT) zx_sl(i) = RLSTT
190    ENDDO
191   
192
193    DO i = 1, knon
194! Q
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)) ) &
201            / zx_oq(i)
202       zx_nq(i) = beta(i) * zx_coefq(i) * (- fqsat * zx_dq_s_dt(i)) &
203            / zx_oq(i)
204       
205! H
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)
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! 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!
260!       if (PRESENT(enth_prec_liq))                   &
261!       enth_prec_liq(i) = rcw * (t1lay(i) - tsurf(i)) * &
262!                          precip_rain(i)
263!       if (PRESENT(enth_prec_sol))                  &
264!       enth_prec_sol(i) = rcs * (t1lay(i) - tsurf(i)) * &
265!                          precip_snow(i)
266! On calcule par rapport a T=0
267       if (PRESENT(sens_prec_liq))                   &
268         sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * &
269                          precip_rain(i)
270       if (PRESENT(sens_prec_sol))                   &
271         sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * &
272                          precip_snow(i)
273       if (PRESENT(lat_prec_liq))                    &
274         lat_prec_liq(i) =  precip_rain(i) * (RLVTT - RLVTT)
275       if (PRESENT(lat_prec_sol))                    &
276         lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
277    ENDDO
278
279   
280!       if (PRESENT(sens_prec_liq))                  &
281!          WRITE(*,*)' calculs_fluxs sens_prec_liq (min, max)', &
282!          MINVAL(sens_prec_liq(1:knon)), MAXVAL(sens_prec_liq(1:knon))
283!       if (PRESENT(sens_prec_sol))                  &
284!          WRITE(*,*)' calculs_fluxs sens_prec_sol (min, max)', &
285!          MINVAL(sens_prec_sol(1:knon)), MAXVAL(sens_prec_sol(1:knon))
286 
287!
288!****************************************************************************************
289!
290  END SUBROUTINE calcul_fluxs
291!
292!****************************************************************************************
293!
294  SUBROUTINE calcul_flux_wind(knon, dtime, &
295       u0, v0, u1, v1, gustiness, cdrag_m, &
296       AcoefU, AcoefV, BcoefU, BcoefV, &
297       p1lay, t1lay, &
298       flux_u1, flux_v1)
299
300    USE dimphy
301    INCLUDE "YOMCST.h"
302    INCLUDE "clesphys.h"
303
304! Input arguments
305!****************************************************************************************
306    INTEGER, INTENT(IN)                  :: knon
307    REAL, INTENT(IN)                     :: dtime
308    REAL, DIMENSION(klon), INTENT(IN)    :: u0, v0  ! u and v at niveau 0
309    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness  ! u and v at niveau 1
310    REAL, DIMENSION(klon), INTENT(IN)    :: cdrag_m ! cdrag pour momentum
311    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
312    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay   ! pression 1er niveau (milieu de couche)
313    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay   ! temperature
314! Output arguments
315!****************************************************************************************
316    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1
317    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v1
318
319! Local variables
320!****************************************************************************************
321    INTEGER                              :: i
322    REAL                                 :: mod_wind, buf
323
324!****************************************************************************************
325! Calculate the surface flux
326!
327!****************************************************************************************
328    DO i=1,knon
329       mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
330       buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
331       flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
332       flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
333    END DO
334
335  END SUBROUTINE calcul_flux_wind
336!
337!****************************************************************************************
338!
339END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.