source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/calcul_fluxs_mod.F90 @ 3363

Last change on this file since 3363 was 3331, checked in by acozic, 6 years ago

Add modification for isotopes

  • Property svn:executable set to *
File size: 13.8 KB
Line 
1!
2! $Id: calcul_fluxs_mod.F90 2538 2016-06-03 14:12:16Z fairhead $
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!        write(*,*) 'calcul_flux 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
134!
135! Initialisation
136!****************************************************************************************
137    evap = 0.
138    fluxsens=0.
139    fluxlat=0.
140    dflux_s = 0.
141    dflux_l = 0.
142    if (PRESENT(sens_prec_liq)) sens_prec_liq = 0.
143    if (PRESENT(sens_prec_sol)) sens_prec_sol = 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).LT.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) .LT. 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!
220! Y'a-t-il fonte de neige?
221!
222!    fonte_neige = (nisurf /= is_oce) .AND. &
223!     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
224!     & .AND. (tsurf_new(i) >= RTT)
225!    if (fonte_neige) tsurf_new(i) = RTT 
226       d_ts(i) = tsurf_new(i) - tsurf(i)
227!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
228!    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
229
230!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
231!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
232       evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
233       fluxlat(i) = - evap(i) * zx_sl(i)
234       fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
235       
236! Derives des flux dF/dTs (W m-2 K-1):
237       dflux_s(i) = zx_nh(i)
238       dflux_l(i) = (zx_sl(i) * zx_nq(i))
239
240! Nouvelle valeure de l'humidite au dessus du sol
241       qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
242       q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
243       qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
244!
245! en cas de fonte de neige
246!
247!    if (fonte_neige) then
248!      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
249!     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
250!     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
251!      bilan_f = max(0., bilan_f)
252!      fq_fonte = bilan_f / zx_sl(i)
253!      snow(i) = max(0., snow(i) - fq_fonte * dtime)
254!      qsol(i) = qsol(i) + (fq_fonte * dtime)
255!    endif
256!!$    if (nisurf == is_ter)  &
257!!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
258!!$    qsol(i) = min(qsol(i), max_eau_sol)
259!
260! calcul de l'enthalpie des precipitations liquides et solides
261!
262!       if (PRESENT(enth_prec_liq))                   &
263!       enth_prec_liq(i) = rcw * (t1lay(i) - tsurf(i)) * &
264!                          precip_rain(i)
265!       if (PRESENT(enth_prec_sol))                  &
266!       enth_prec_sol(i) = rcs * (t1lay(i) - tsurf(i)) * &
267!                          precip_snow(i)
268! On calcule par rapport a T=0
269       if (PRESENT(sens_prec_liq))                   &
270         sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * &
271                          precip_rain(i)
272       if (PRESENT(sens_prec_sol))                   &
273         sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * &
274                          precip_snow(i)
275       if (PRESENT(lat_prec_liq))                    &
276         lat_prec_liq(i) =  precip_rain(i) * (RLVTT - RLVTT)
277       if (PRESENT(lat_prec_sol))                    &
278         lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
279    ENDDO
280
281   
282!       if (PRESENT(sens_prec_liq))                  &
283!          WRITE(*,*)' calculs_fluxs sens_prec_liq (min, max)', &
284!          MINVAL(sens_prec_liq(1:knon)), MAXVAL(sens_prec_liq(1:knon))
285!       if (PRESENT(sens_prec_sol))                  &
286!          WRITE(*,*)' calculs_fluxs sens_prec_sol (min, max)', &
287!          MINVAL(sens_prec_sol(1:knon)), MAXVAL(sens_prec_sol(1:knon))
288 
289!
290!****************************************************************************************
291!
292  END SUBROUTINE calcul_fluxs
293!
294!****************************************************************************************
295!
296  SUBROUTINE calcul_flux_wind(knon, dtime, &
297       u0, v0, u1, v1, gustiness, cdrag_m, &
298       AcoefU, AcoefV, BcoefU, BcoefV, &
299       p1lay, t1lay, &
300       flux_u1, flux_v1)
301
302    USE dimphy
303    INCLUDE "YOMCST.h"
304    INCLUDE "clesphys.h"
305
306! Input arguments
307!****************************************************************************************
308    INTEGER, INTENT(IN)                  :: knon
309    REAL, INTENT(IN)                     :: dtime
310    REAL, DIMENSION(klon), INTENT(IN)    :: u0, v0  ! u and v at niveau 0
311    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness  ! u and v at niveau 1
312    REAL, DIMENSION(klon), INTENT(IN)    :: cdrag_m ! cdrag pour momentum
313    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
314    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay   ! pression 1er niveau (milieu de couche)
315    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay   ! temperature
316! Output arguments
317!****************************************************************************************
318    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1
319    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v1
320
321! Local variables
322!****************************************************************************************
323    INTEGER                              :: i
324    REAL                                 :: mod_wind, buf
325
326!****************************************************************************************
327! Calculate the surface flux
328!
329!****************************************************************************************
330    DO i=1,knon
331       mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
332       buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
333       flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
334       flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
335    END DO
336
337  END SUBROUTINE calcul_flux_wind
338!
339!****************************************************************************************
340!
341END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.