/[lmdze]/trunk/phylmd/clqh.f
ViewVC logotype

Annotation of /trunk/phylmd/clqh.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 1 month ago) by guez
File size: 11607 byte(s)
Changed all ".f90" suffixes to ".f".
1 guez 49 module clqh_m
2 guez 3
3 guez 38 IMPLICIT none
4 guez 3
5 guez 49 contains
6 guez 3
7 guez 72 SUBROUTINE clqh(dtime, itime, jour, debut, rlat, knon, nisurf, knindex, &
8     pctsrf, soil_model, tsoil, qsol, ok_veget, ocean, rmu0, co2_ppm, &
9     rugos, rugoro, u1lay, v1lay, coef, t, q, ts, paprs, pplay, delp, &
10     radsol, albedo, alblw, snow, qsurf, precip_rain, precip_snow, fder, &
11     swnet, fluxlat, pctsrf_new, agesno, d_t, d_q, d_ts, z0_new, flux_t, &
12     flux_q, dflux_s, dflux_l, fqcalving, ffonte, run_off_lic_0, flux_o, &
13     flux_g, tslab, seaice)
14 guez 3
15 guez 62 ! Author: Z. X. Li (LMD/CNRS)
16 guez 49 ! Date: 1993/08/18
17     ! Objet : diffusion verticale de "q" et de "h"
18 guez 3
19 guez 49 USE conf_phys_m, ONLY : iflag_pbl
20     USE dimens_m, ONLY : iim, jjm
21     USE dimphy, ONLY : klev, klon, zmasq
22     USE dimsoil, ONLY : nsoilmx
23     USE indicesol, ONLY : is_ter, nbsrf
24 guez 54 USE interfsurf_hq_m, ONLY : interfsurf_hq
25 guez 49 USE suphec_m, ONLY : rcpd, rd, rg, rkappa
26 guez 38
27 guez 49 ! Arguments:
28 guez 70 INTEGER, intent(in):: knon
29 guez 49 REAL, intent(in):: dtime ! intervalle du temps (s)
30     REAL u1lay(klon) ! vitesse u de la 1ere couche (m/s)
31     REAL v1lay(klon) ! vitesse v de la 1ere couche (m/s)
32 guez 70
33     REAL, intent(in):: coef(:, :) ! (knon, klev)
34     ! Le coefficient d'echange (m**2/s) multiplie par le cisaillement
35     ! du vent (dV/dz). La premiere valeur indique la valeur de Cdrag
36     ! (sans unite).
37    
38 guez 49 REAL t(klon, klev) ! temperature (K)
39     REAL q(klon, klev) ! humidite specifique (kg/kg)
40     REAL ts(klon) ! temperature du sol (K)
41     REAL evap(klon) ! evaporation au sol
42     REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
43     REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
44     REAL delp(klon, klev) ! epaisseur de couche en pression (Pa)
45     REAL radsol(klon) ! ray. net au sol (Solaire+IR) W/m2
46     REAL albedo(klon) ! albedo de la surface
47     REAL alblw(klon)
48     REAL snow(klon) ! hauteur de neige
49     REAL qsurf(klon) ! humidite de l'air au dessus de la surface
50     real precip_rain(klon), precip_snow(klon)
51     REAL agesno(klon)
52     REAL rugoro(klon)
53     REAL run_off_lic_0(klon)! runof glacier au pas de temps precedent
54 guez 62 integer, intent(in):: jour ! jour de l'annee en cours
55     real, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal
56 guez 49 real rugos(klon) ! rugosite
57     integer knindex(klon)
58 guez 62 real, intent(in):: pctsrf(klon, nbsrf)
59 guez 72 real, intent(in):: rlat(klon)
60 guez 49 logical ok_veget
61     REAL co2_ppm ! taux CO2 atmosphere
62     character(len=*), intent(in):: ocean
63 guez 3
64 guez 49 REAL d_t(klon, klev) ! incrementation de "t"
65     REAL d_q(klon, klev) ! incrementation de "q"
66     REAL d_ts(klon) ! incrementation de "ts"
67     REAL flux_t(klon, klev) ! (diagnostic) flux de la chaleur
68     ! sensible, flux de Cp*T, positif vers
69     ! le bas: j/(m**2 s) c.a.d.: W/m2
70     REAL flux_q(klon, klev) ! flux de la vapeur d'eau:kg/(m**2 s)
71     REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
72     REAL dflux_l(klon) ! derivee du flux latent dF/dTs
73     !IM cf JLD
74     ! Flux thermique utiliser pour fondre la neige
75     REAL ffonte(klon)
76     ! Flux d'eau "perdue" par la surface et nĂ©cessaire pour que limiter la
77     ! hauteur de neige, en kg/m2/s
78     REAL fqcalving(klon)
79     !IM "slab" ocean
80     REAL tslab(klon) !temperature du slab ocean (K) (OCEAN='slab ')
81     REAL seaice(klon) ! glace de mer en kg/m2
82     REAL flux_o(klon) ! flux entre l'ocean et l'atmosphere W/m2
83     REAL flux_g(klon) ! flux entre l'ocean et la glace de mer W/m2
84 guez 3
85 guez 49 INTEGER i, k
86     REAL zx_cq(klon, klev)
87     REAL zx_dq(klon, klev)
88     REAL zx_ch(klon, klev)
89     REAL zx_dh(klon, klev)
90     REAL zx_buf1(klon)
91     REAL zx_buf2(klon)
92     REAL zx_coef(klon, klev)
93     REAL local_h(klon, klev) ! enthalpie potentielle
94     REAL local_q(klon, klev)
95     REAL local_ts(klon)
96     REAL psref(klon) ! pression de reference pour temperature potent.
97     REAL zx_pkh(klon, klev), zx_pkf(klon, klev)
98 guez 3
99 guez 49 ! contre-gradient pour la vapeur d'eau: (kg/kg)/metre
100     REAL gamq(klon, 2:klev)
101     ! contre-gradient pour la chaleur sensible: Kelvin/metre
102     REAL gamt(klon, 2:klev)
103     REAL z_gamaq(klon, 2:klev), z_gamah(klon, 2:klev)
104     REAL zdelz
105 guez 3
106 guez 49 ! Rajout pour l'interface
107     integer, intent(in):: itime
108     integer nisurf
109     logical, intent(in):: debut
110     real zlev1(klon)
111 guez 72 real fder(klon)
112 guez 49 real temp_air(klon), spechum(klon)
113     real epot_air(klon), ccanopy(klon)
114     real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
115     real petBcoef(klon), peqBcoef(klon)
116 guez 72 real swnet(klon), swdown(klon)
117 guez 49 real p1lay(klon)
118     !$$$C PB ajout pour soil
119     LOGICAL, intent(in):: soil_model
120     REAL tsoil(klon, nsoilmx)
121     REAL qsol(klon)
122 guez 3
123 guez 49 ! Parametres de sortie
124     real fluxsens(klon), fluxlat(klon)
125 guez 72 real tsurf_new(klon), alb_new(klon)
126     real z0_new(klon)
127 guez 49 real pctsrf_new(klon, nbsrf)
128     ! JLD
129     real zzpk
130 guez 3
131 guez 49 character (len = 20) :: modname = 'Debut clqh'
132     LOGICAL check
133     PARAMETER (check=.false.)
134 guez 38
135 guez 49 !----------------------------------------------------------------
136 guez 3
137 guez 49 if (check) THEN
138     write(*, *) modname, ' nisurf=', nisurf
139     !C call flush(6)
140     endif
141    
142     if (check) THEN
143     WRITE(*, *)' qsurf (min, max)' &
144     , minval(qsurf(1:knon)), maxval(qsurf(1:knon))
145     !C call flush(6)
146     ENDIF
147    
148     if (iflag_pbl.eq.1) then
149     do k = 3, klev
150     do i = 1, knon
151     gamq(i, k)= 0.0
152     gamt(i, k)= -1.0e-03
153     enddo
154     enddo
155     do i = 1, knon
156     gamq(i, 2) = 0.0
157     gamt(i, 2) = -2.5e-03
158     enddo
159     else
160     do k = 2, klev
161     do i = 1, knon
162     gamq(i, k) = 0.0
163     gamt(i, k) = 0.0
164     enddo
165     enddo
166     endif
167    
168     DO i = 1, knon
169     psref(i) = paprs(i, 1) !pression de reference est celle au sol
170     local_ts(i) = ts(i)
171     ENDDO
172     DO k = 1, klev
173     DO i = 1, knon
174     zx_pkh(i, k) = (psref(i)/paprs(i, k))**RKAPPA
175     zx_pkf(i, k) = (psref(i)/pplay(i, k))**RKAPPA
176     local_h(i, k) = RCPD * t(i, k) * zx_pkf(i, k)
177     local_q(i, k) = q(i, k)
178     ENDDO
179     ENDDO
180    
181     ! Convertir les coefficients en variables convenables au calcul:
182    
183     DO k = 2, klev
184     DO i = 1, knon
185     zx_coef(i, k) = coef(i, k)*RG/(pplay(i, k-1)-pplay(i, k)) &
186     *(paprs(i, k)*2/(t(i, k)+t(i, k-1))/RD)**2
187     zx_coef(i, k) = zx_coef(i, k) * dtime*RG
188     ENDDO
189     ENDDO
190    
191     ! Preparer les flux lies aux contre-gardients
192    
193     DO k = 2, klev
194     DO i = 1, knon
195     zdelz = RD * (t(i, k-1)+t(i, k))/2.0 / RG /paprs(i, k) &
196     *(pplay(i, k-1)-pplay(i, k))
197     z_gamaq(i, k) = gamq(i, k) * zdelz
198     z_gamah(i, k) = gamt(i, k) * zdelz *RCPD * zx_pkh(i, k)
199     ENDDO
200     ENDDO
201     DO i = 1, knon
202     zx_buf1(i) = zx_coef(i, klev) + delp(i, klev)
203     zx_cq(i, klev) = (local_q(i, klev)*delp(i, klev) &
204     -zx_coef(i, klev)*z_gamaq(i, klev))/zx_buf1(i)
205     zx_dq(i, klev) = zx_coef(i, klev) / zx_buf1(i)
206    
207     zzpk=(pplay(i, klev)/psref(i))**RKAPPA
208     zx_buf2(i) = zzpk*delp(i, klev) + zx_coef(i, klev)
209     zx_ch(i, klev) = (local_h(i, klev)*zzpk*delp(i, klev) &
210     -zx_coef(i, klev)*z_gamah(i, klev))/zx_buf2(i)
211     zx_dh(i, klev) = zx_coef(i, klev) / zx_buf2(i)
212     ENDDO
213     DO k = klev-1, 2 , -1
214     DO i = 1, knon
215     zx_buf1(i) = delp(i, k)+zx_coef(i, k) &
216     +zx_coef(i, k+1)*(1.-zx_dq(i, k+1))
217     zx_cq(i, k) = (local_q(i, k)*delp(i, k) &
218     +zx_coef(i, k+1)*zx_cq(i, k+1) &
219     +zx_coef(i, k+1)*z_gamaq(i, k+1) &
220     -zx_coef(i, k)*z_gamaq(i, k))/zx_buf1(i)
221     zx_dq(i, k) = zx_coef(i, k) / zx_buf1(i)
222    
223     zzpk=(pplay(i, k)/psref(i))**RKAPPA
224     zx_buf2(i) = zzpk*delp(i, k)+zx_coef(i, k) &
225     +zx_coef(i, k+1)*(1.-zx_dh(i, k+1))
226     zx_ch(i, k) = (local_h(i, k)*zzpk*delp(i, k) &
227     +zx_coef(i, k+1)*zx_ch(i, k+1) &
228     +zx_coef(i, k+1)*z_gamah(i, k+1) &
229     -zx_coef(i, k)*z_gamah(i, k))/zx_buf2(i)
230     zx_dh(i, k) = zx_coef(i, k) / zx_buf2(i)
231     ENDDO
232     ENDDO
233    
234     DO i = 1, knon
235     zx_buf1(i) = delp(i, 1) + zx_coef(i, 2)*(1.-zx_dq(i, 2))
236     zx_cq(i, 1) = (local_q(i, 1)*delp(i, 1) &
237     +zx_coef(i, 2)*(z_gamaq(i, 2)+zx_cq(i, 2))) &
238     /zx_buf1(i)
239     zx_dq(i, 1) = -1. * RG / zx_buf1(i)
240    
241     zzpk=(pplay(i, 1)/psref(i))**RKAPPA
242     zx_buf2(i) = zzpk*delp(i, 1) + zx_coef(i, 2)*(1.-zx_dh(i, 2))
243     zx_ch(i, 1) = (local_h(i, 1)*zzpk*delp(i, 1) &
244     +zx_coef(i, 2)*(z_gamah(i, 2)+zx_ch(i, 2))) &
245     /zx_buf2(i)
246     zx_dh(i, 1) = -1. * RG / zx_buf2(i)
247     ENDDO
248    
249     ! Appel a interfsurf (appel generique) routine d'interface avec la surface
250    
251     ! initialisation
252     petAcoef =0.
253     peqAcoef = 0.
254     petBcoef =0.
255     peqBcoef = 0.
256     p1lay =0.
257    
258     petAcoef(1:knon) = zx_ch(1:knon, 1)
259     peqAcoef(1:knon) = zx_cq(1:knon, 1)
260     petBcoef(1:knon) = zx_dh(1:knon, 1)
261     peqBcoef(1:knon) = zx_dq(1:knon, 1)
262 guez 70 tq_cdrag(1:knon) =coef(:knon, 1)
263 guez 49 temp_air(1:knon) =t(1:knon, 1)
264     epot_air(1:knon) =local_h(1:knon, 1)
265     spechum(1:knon)=q(1:knon, 1)
266     p1lay(1:knon) = pplay(1:knon, 1)
267     zlev1(1:knon) = delp(1:knon, 1)
268    
269     if(nisurf.eq.is_ter) THEN
270     swdown(1:knon) = swnet(1:knon)/(1-albedo(1:knon))
271     else
272     swdown(1:knon) = swnet(1:knon)
273     endif
274     ccanopy = co2_ppm
275    
276 guez 72 CALL interfsurf_hq(itime, dtime, jour, rmu0, klon, iim, jjm, &
277     nisurf, knon, knindex, pctsrf, rlat, debut, &
278     ok_veget, soil_model, nsoilmx, tsoil, qsol, u1lay, v1lay, &
279     temp_air, spechum, tq_cdrag, petAcoef, peqAcoef, &
280     petBcoef, peqBcoef, precip_rain, precip_snow, &
281     fder, rugos, rugoro, &
282     snow, qsurf, ts, p1lay, psref, radsol, ocean, &
283     evap, fluxsens, fluxlat, dflux_l, dflux_s, tsurf_new, &
284     alb_new, alblw, z0_new, pctsrf_new, agesno, fqcalving, &
285     ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)
286 guez 49
287     do i = 1, knon
288     flux_t(i, 1) = fluxsens(i)
289     flux_q(i, 1) = - evap(i)
290     d_ts(i) = tsurf_new(i) - ts(i)
291     albedo(i) = alb_new(i)
292     enddo
293    
294     !==== une fois on a zx_h_ts, on peut faire l'iteration ========
295     DO i = 1, knon
296     local_h(i, 1) = zx_ch(i, 1) + zx_dh(i, 1)*flux_t(i, 1)*dtime
297     local_q(i, 1) = zx_cq(i, 1) + zx_dq(i, 1)*flux_q(i, 1)*dtime
298     ENDDO
299     DO k = 2, klev
300     DO i = 1, knon
301     local_q(i, k) = zx_cq(i, k) + zx_dq(i, k)*local_q(i, k-1)
302     local_h(i, k) = zx_ch(i, k) + zx_dh(i, k)*local_h(i, k-1)
303     ENDDO
304     ENDDO
305     !======================================================================
306     !== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
307     !== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
308     DO k = 2, klev
309     DO i = 1, knon
310     flux_q(i, k) = (zx_coef(i, k)/RG/dtime) &
311     * (local_q(i, k)-local_q(i, k-1)+z_gamaq(i, k))
312     flux_t(i, k) = (zx_coef(i, k)/RG/dtime) &
313     * (local_h(i, k)-local_h(i, k-1)+z_gamah(i, k)) &
314     / zx_pkh(i, k)
315     ENDDO
316     ENDDO
317     !======================================================================
318     ! Calcul tendances
319     DO k = 1, klev
320     DO i = 1, knon
321     d_t(i, k) = local_h(i, k)/zx_pkf(i, k)/RCPD - t(i, k)
322     d_q(i, k) = local_q(i, k) - q(i, k)
323     ENDDO
324     ENDDO
325    
326     END SUBROUTINE clqh
327    
328     end module clqh_m

  ViewVC Help
Powered by ViewVC 1.1.21