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

Annotation of /trunk/phylmd/clqh.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 12 - (hide annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/clqh.f
File size: 13454 byte(s)
-- Minor modification of input/output:

Created procedure "read_logic". Variables of module "logic" are read
by "read_logic" instead of "conf_gcm". Variable "offline" of module
"conf_gcm" is read from namelist instead of "*.def".

Deleted arguments "dtime", "co2_ppm_etat0", "solaire_etat0",
"tabcntr0" and local variables "radpas", "tab_cntrl" of
"phyetat0". "phyetat0" does not read "controle" in "startphy.nc" any
longer. "phyetat0" now reads global attribute "itau_phy" from
"startphy.nc". "phyredem" does not create variable "controle" in
"startphy.nc" any longer. "phyredem" now writes global attribute
"itau_phy" of "startphy.nc". Deleted argument "tabcntr0" of
"printflag". Removed diagnostic messages written by "printflag" for
comparison of the variable "controle" of "startphy.nc" and the
variables read from "*.def" or namelist input.

-- Removing unwanted functionality:

Removed variable "lunout" from module "iniprint", replaced everywhere
by standard output.

Removed case "ocean == 'couple'" in "clmain", "interfsurf_hq" and
"physiq". Removed procedure "interfoce_cpl".

-- Should not change anything at run time:

Automated creation of graphs in documentation. More documentation on
input files.

Converted Fortran files to free format: "phyredem.f90", "printflag.f90".

Split module "clesphy" into "clesphys" and "clesphys2".

Removed variables "conser", "leapf", "forward", "apphys", "apdiss" and
"statcl" from module "logic". Added arguments "conser" to "advect",
"leapf" to "integrd". Added local variables "forward", "leapf",
"apphys", "conser", "apdiss" in "leapfrog".

Added intent attributes.

Deleted arguments "dtime" of "phyredem", "pdtime" of "flxdtdq", "sh"
of "phytrac", "dt" of "yamada".

Deleted local variables "dtime", "co2_ppm_etat0", "solaire_etat0",
"length", "tabcntr0" in "physiq". Replaced all references to "dtime"
by references to "pdtphys".

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 guez 12 REAL, intent(in):: dtime ! intervalle du temps (s)
42 guez 3 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 guez 12 character(len=*), intent(in):: ocean
75 guez 3 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 guez 7 integer, intent(in):: itime
135 guez 3 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 guez 12 LOGICAL, intent(in):: soil_model
148 guez 3 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