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

Annotation of /trunk/phylmd/clqh.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21