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

Annotation of /trunk/Sources/phylmd/clqh.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 62 - (hide annotations)
Thu Jul 26 14:37:37 2012 UTC (11 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/clqh.f90
File size: 12671 byte(s)
Changed handling of compiler in compilation system.

Removed the prefix letters "y", "p", "t" or "z" in some names of variables.

Replaced calls to NetCDF by calls to NetCDF95.

Extracted "ioget_calendar" procedures from "calendar.f90" into a
separate file.

Extracted to a separate file, "mathop2.f90", procedures that were not
part of the generic interface "mathop" in "mathop.f90".

Removed computation of "dq" in "bilan_dyn", which was not used.

In "iniadvtrac", removed schemes 20 Slopes and 30 Prather. Was not
compatible with declarations of array sizes.

In "clcdrag", "ustarhb", "vdif_kcay", "yamada4" and "coefkz", changed
the size of some arrays from "klon" to "knon".

Removed possible call to "conema3" in "physiq".

Removed unused argument "cd" in "yamada".

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

  ViewVC Help
Powered by ViewVC 1.1.21