source: LMDZ4/trunk/libf/phy_IPCC_AR4/calcul_fluxs_mod.F90 @ 945

Last change on this file since 945 was 868, checked in by Laurent Fairhead, 17 years ago

Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4
par la version de la physique avec thermique. On garde le repertoire phylmd
pour un petit moment pour que les utilisateurs ne soient pas trop perdus ...
phy_IPCC_AR4 = phylmd
LF
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 9.0 KB
Line 
1!
2! $Header$
3!
4MODULE calcul_fluxs_mod
5
6
7CONTAINS
8  SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
9       tsurf, p1lay, cal, beta, coef1lay, ps, &
10       precip_rain, precip_snow, snow, qsurf, &
11       radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
12       petAcoef, peqAcoef, petBcoef, peqBcoef, &
13       tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
14   
15    USE dimphy, ONLY : klon
16
17! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
18! une temperature de surface (au cas ou ok_veget = false)
19!
20! L. Fairhead 4/2000
21!
22! input:
23!   knon         nombre de points a traiter
24!   nisurf       surface a traiter
25!   tsurf        temperature de surface
26!   p1lay        pression 1er niveau (milieu de couche)
27!   cal          capacite calorifique du sol
28!   beta         evap reelle
29!   coef1lay     coefficient d'echange
30!   ps           pression au sol
31!   precip_rain  precipitations liquides
32!   precip_snow  precipitations solides
33!   snow         champs hauteur de neige
34!   runoff       runoff en cas de trop plein
35!   petAcoef     coeff. A de la resolution de la CL pour t
36!   peqAcoef     coeff. A de la resolution de la CL pour q
37!   petBcoef     coeff. B de la resolution de la CL pour t
38!   peqBcoef     coeff. B de la resolution de la CL pour q
39!   radsol       rayonnement net aus sol (LW + SW)
40!   dif_grnd     coeff. diffusion vers le sol profond
41!
42! output:
43!   tsurf_new    temperature au sol
44!   qsurf        humidite de l'air au dessus du sol
45!   fluxsens     flux de chaleur sensible
46!   fluxlat      flux de chaleur latente
47!   dflux_s      derivee du flux de chaleur sensible / Ts
48!   dflux_l      derivee du flux de chaleur latente  / Ts
49!
50
51    INCLUDE "YOETHF.h"
52    INCLUDE "FCTTRE.h"
53    INCLUDE "indicesol.h"
54    INCLUDE "YOMCST.h"
55
56! Parametres d'entree
57!****************************************************************************************
58    INTEGER, INTENT(IN)                  :: knon, nisurf
59    REAL   , INTENT(IN)                  :: dtime
60    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
61    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
62    REAL, DIMENSION(klon), INTENT(IN)    :: ps, q1lay
63    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf, p1lay, cal, beta, coef1lay
64    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow ! pas utiles
65    REAL, DIMENSION(klon), INTENT(IN)    :: radsol, dif_grnd
66    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay, u1lay, v1lay
67
68! Parametres entree-sorties
69!****************************************************************************************
70    REAL, DIMENSION(klon), INTENT(INOUT) :: snow  ! snow pas utile
71
72! Parametres sorties
73!****************************************************************************************
74    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
75    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new, evap, fluxsens, fluxlat
76    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l
77
78! Variables locales
79!****************************************************************************************
80    INTEGER                              :: i
81    REAL, DIMENSION(klon)                :: zx_mh, zx_nh, zx_oh
82    REAL, DIMENSION(klon)                :: zx_mq, zx_nq, zx_oq
83    REAL, DIMENSION(klon)                :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
84    REAL, DIMENSION(klon)                :: zx_sl, zx_k1
85    REAL, DIMENSION(klon)                :: d_ts
86    REAL                                 :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
87    REAL                                 :: qsat_new, q1_new
88    REAL, PARAMETER                      :: t_grnd = 271.35, t_coup = 273.15
89    REAL, PARAMETER                      :: max_eau_sol = 150.0
90    CHARACTER (len = 20)                 :: modname = 'calcul_fluxs'
91    LOGICAL                              :: fonte_neige
92    LOGICAL, SAVE                        :: check = .FALSE.
93    !$OMP THREADPRIVATE(check)
94
95! End definition
96!****************************************************************************************
97
98    IF (check) WRITE(*,*)'Entree ', modname,' surface = ',nisurf
99   
100    IF (check) THEN
101       WRITE(*,*)' radsol (min, max)', &
102            MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
103       CALL flush(6)
104    ENDIF
105 
106! Traitement neige et humidite du sol
107!****************************************************************************************
108!
109!!$  WRITE(*,*)'test calcul_flux, surface ', nisurf
110!!PB test
111!!$    if (nisurf == is_oce) then
112!!$      snow = 0.
113!!$      qsol = max_eau_sol
114!!$    else
115!!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
116!!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
117!!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)
118!!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
119!!$    endif
120!!$    IF (nisurf /= is_ter) qsol = max_eau_sol
121
122
123!
124! Initialisation
125!****************************************************************************************
126    evap = 0.
127    fluxsens=0.
128    fluxlat=0.
129    dflux_s = 0.
130    dflux_l = 0.       
131!
132! zx_qs = qsat en kg/kg
133!****************************************************************************************
134    DO i = 1, knon
135       zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
136       IF (thermcep) THEN
137          zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
138          zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
139          zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
140          zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
141          zx_qs=MIN(0.5,zx_qs)
142          zcor=1./(1.-retv*zx_qs)
143          zx_qs=zx_qs*zcor
144          zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
145               /RLVTT / zx_pkh(i)
146       ELSE
147          IF (tsurf(i).LT.t_coup) THEN
148             zx_qs = qsats(tsurf(i)) / ps(i)
149             zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
150                  / zx_pkh(i)
151          ELSE
152             zx_qs = qsatl(tsurf(i)) / ps(i)
153             zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
154                  / zx_pkh(i)
155          ENDIF
156       ENDIF
157       zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
158       zx_qsat(i) = zx_qs
159       zx_coef(i) = coef1lay(i) * &
160            (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) * &
161            p1lay(i)/(RD*t1lay(i))
162       
163    ENDDO
164
165
166! === Calcul de la temperature de surface ===
167! zx_sl = chaleur latente d'evaporation ou de sublimation
168!****************************************************************************************
169
170    DO i = 1, knon
171       zx_sl(i) = RLVTT
172       IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
173       zx_k1(i) = zx_coef(i)
174    ENDDO
175   
176
177    DO i = 1, knon
178! Q
179       zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
180       zx_mq(i) = beta(i) * zx_k1(i) * &
181            (peqAcoef(i) - zx_qsat(i) + &
182            zx_dq_s_dt(i) * tsurf(i)) &
183            / zx_oq(i)
184       zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
185            / zx_oq(i)
186       
187! H
188       zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
189       zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
190       zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
191     
192! Tsurface
193       tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
194            (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
195            + dif_grnd(i) * t_grnd * dtime)/ &
196            ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
197            zx_nh(i) + zx_sl(i) * zx_nq(i)) & 
198            + dtime * dif_grnd(i))
199
200!
201! Y'a-t-il fonte de neige?
202!
203!    fonte_neige = (nisurf /= is_oce) .AND. &
204!     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
205!     & .AND. (tsurf_new(i) >= RTT)
206!    if (fonte_neige) tsurf_new(i) = RTT 
207       d_ts(i) = tsurf_new(i) - tsurf(i)
208!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
209!    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
210
211!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
212!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
213       evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
214       fluxlat(i) = - evap(i) * zx_sl(i)
215       fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
216       
217! Derives des flux dF/dTs (W m-2 K-1):
218       dflux_s(i) = zx_nh(i)
219       dflux_l(i) = (zx_sl(i) * zx_nq(i))
220
221! Nouvelle valeure de l'humidite au dessus du sol
222       qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
223       q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
224       qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
225!
226! en cas de fonte de neige
227!
228!    if (fonte_neige) then
229!      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
230!     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
231!     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
232!      bilan_f = max(0., bilan_f)
233!      fq_fonte = bilan_f / zx_sl(i)
234!      snow(i) = max(0., snow(i) - fq_fonte * dtime)
235!      qsol(i) = qsol(i) + (fq_fonte * dtime)
236!    endif
237!!$    if (nisurf == is_ter)  &
238!!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
239!!$    qsol(i) = min(qsol(i), max_eau_sol)
240    ENDDO
241!
242!****************************************************************************************
243!
244  END SUBROUTINE calcul_fluxs
245
246END MODULE calcul_fluxs_mod
Note: See TracBrowser for help on using the repository browser.