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

Annotation of /trunk/phylmd/clqh.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 70 - (hide annotations)
Mon Jun 24 15:39:52 2013 UTC (10 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/clqh.f90
File size: 12389 byte(s)
In procedure, "addfi" access directly the module variable "dtphys"
instead of going through an argument.

In "conflx", do not create a local variable for temperature with
reversed order of vertical levels. Instead, give an actual argument
with reversed order in "physiq".

Changed names of variables "rmd" and "rmv" from module "suphec_m" to
"md" and "mv".

In "hgardfou", print only the first temperature out of range found.

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

  ViewVC Help
Powered by ViewVC 1.1.21