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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 54 - (hide annotations)
Tue Dec 6 15:07:04 2011 UTC (12 years, 4 months ago) by guez
Original Path: trunk/libf/phylmd/clqh.f90
File size: 12615 byte(s)
Removed Numerical Recipes procedure "ran1". Replaced calls to "ran1"
in "inidissip" by calls to intrinsic procedures.

Split file "interface_surf.f90" into a file with a module containing
only variables, "interface_surf", and single-procedure files. Gathered
files into directory "Interface_surf".

Added argument "cdivu" to "gradiv" and "gradiv2", "cdivh" to
"divgrad2" and "divgrad", and "crot" to "nxgraro2" and
"nxgrarot". "dissip" now uses variables "cdivu", "cdivh" and "crot"
from module "inidissip_m", so it can pass them to "gradiv2",
etc. Thanks to this modification, we avoid a circular dependency
betwwen "inidissip.f90" and "gradiv2.f90", etc. The value -1. used by
"gradiv2", for instance, during computation of eigenvalues is not the
value "cdivu" computed by "inidissip".

Extracted procedure "start_inter_3d" from module "startdyn", to its
own module.

In "inidissip", unrolled loop on "ii". I find it clearer now.

Moved variables "matriceun", "matriceus", "matricevn", "matricevs",
"matrinvn" and "matrinvs" from module "parafilt" to module
"inifilr_m". Moved variables "jfiltnu", "jfiltnv", "jfiltsu",
"jfiltsv" from module "coefils" to module "inifilr_m".

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 49 ! Author: Z.X. Li (LMD/CNRS)!
17     ! 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     real date0
32     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     integer jour ! jour de l'annee en cours
55     real rmu0(klon) ! cosinus de l'angle solaire zenithal
56     real rugos(klon) ! rugosite
57     integer knindex(klon)
58     real pctsrf(klon, nbsrf)
59     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