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

Annotation of /trunk/phylmd/clqh.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 38 - (hide annotations)
Thu Jan 6 17:52:19 2011 UTC (13 years, 3 months ago) by guez
Original Path: trunk/libf/phylmd/clqh.f90
File size: 12235 byte(s)
Extracted ASCII art from "inigeom" into a separate text file in the
documentation.

"test_disvert" now creates a separate file for layer thicknesses.

Moved variables from module "yomcst" to module "suphec_m" because this
is where those variables are defined. Kept in "yomcst" only parameters
of Earth orbit. Gave the attribute "parameter" to some variables of
module "suphec_m".

Variables of module "yoethf" were defined in procedure "suphec". Moved
these definitions to a new procedure "yoethf" in module "yoethf_m".

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

  ViewVC Help
Powered by ViewVC 1.1.21