/[lmdze]/trunk/libf/phylmd/diagphy.f
ViewVC logotype

Annotation of /trunk/libf/phylmd/diagphy.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 10 - (hide annotations)
Fri Apr 18 14:45:53 2008 UTC (16 years, 1 month ago) by guez
File size: 14253 byte(s)
Added NetCDF directory "/home/guez/include" in "g95.mk" and
"nag_tools.mk".

Added some "intent" attributes in "PVtheta", "advtrac", "caladvtrac",
"calfis", "diagedyn", "dissip", "vlspltqs", "aeropt", "ajsec",
"calltherm", "clmain", "cltrac", "cltracrn", "concvl", "conema3",
"conflx", "fisrtilp", "newmicro", "nuage", "diagcld1", "diagcld2",
"drag_noro", "lift_noro", "SUGWD", "physiq", "phytrac", "radlwsw", "thermcell".

Removed the case "ierr == 0" in "abort_gcm"; moved call to "histclo"
and messages for end of run from "abort_gcm" to "gcm"; replaced call
to "abort_gcm" in "leapfrog" by exit from outer loop.

In "calfis": removed argument "pp" and variable "unskap"; changed
"pksurcp" from scalar to rank 2; use "pressure_var"; rewrote
computation of "zplev", "zplay", "ztfi", "pcvgt" using "dyn_phy";
added computation of "pls".

Removed unused variable in "dynredem0".

In "exner_hyb": changed "dellta" from scalar to rank 1; replaced call
to "ssum" by call to "sum"; removed variables "xpn" and "xps";
replaced some loops by array expressions.

In "leapfrog": use "pressure_var"; deleted variables "p", "longcles".

Converted common blocks "YOECUMF", "YOEGWD" to modules.

Removed argument "pplay" in "cvltr", "diagetpq", "nflxtr".

Created module "raddimlw" from include file "raddimlw.h".

Corrected call to "new_unit" in "test_disvert".

1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/diagphy.F,v 1.1.1.1 2004/05/19 12:53:08 lmdzadmin Exp $
3     !
4     SUBROUTINE diagphy(airephy,tit,iprt
5     $ , tops, topl, sols, soll, sens
6     $ , evap, rain_fall, snow_fall, ts
7     $ , d_etp_tot, d_qt_tot, d_ec_tot
8     $ , fs_bound, fq_bound)
9     C======================================================================
10     C
11     C Purpose:
12     C Compute the thermal flux and the watter mass flux at the atmosphere
13     c boundaries. Print them and also the atmospheric enthalpy change and
14     C the atmospheric mass change.
15     C
16     C Arguments:
17     C airephy-------input-R- grid area
18     C tit---------input-A15- Comment to be added in PRINT (CHARACTER*15)
19     C iprt--------input-I- PRINT level ( <=0 : no PRINT)
20     C tops(klon)--input-R- SW rad. at TOA (W/m2), positive up.
21     C topl(klon)--input-R- LW rad. at TOA (W/m2), positive down
22     C sols(klon)--input-R- Net SW flux above surface (W/m2), positive up
23     C (i.e. -1 * flux absorbed by the surface)
24     C soll(klon)--input-R- Net LW flux above surface (W/m2), positive up
25     C (i.e. flux emited - flux absorbed by the surface)
26     C sens(klon)--input-R- Sensible Flux at surface (W/m2), positive down
27     C evap(klon)--input-R- Evaporation + sublimation watter vapour mass flux
28     C (kg/m2/s), positive up
29     C rain_fall(klon)
30     C --input-R- Liquid watter mass flux (kg/m2/s), positive down
31     C snow_fall(klon)
32     C --input-R- Solid watter mass flux (kg/m2/s), positive down
33     C ts(klon)----input-R- Surface temperature (K)
34     C d_etp_tot---input-R- Heat flux equivalent to atmospheric enthalpy
35     C change (W/m2)
36     C d_qt_tot----input-R- Mass flux equivalent to atmospheric watter mass
37     C change (kg/m2/s)
38     C d_ec_tot----input-R- Flux equivalent to atmospheric cinetic energy
39     C change (W/m2)
40     C
41     C fs_bound---output-R- Thermal flux at the atmosphere boundaries (W/m2)
42     C fq_bound---output-R- Watter mass flux at the atmosphere boundaries (kg/m2/s)
43     C
44     C J.L. Dufresne, July 2002
45     C======================================================================
46     C
47     use dimens_m
48     use dimphy
49     use YOMCST
50     use yoethf
51     implicit none
52    
53     C
54     C Input variables
55     real airephy(klon)
56     CHARACTER*15 tit
57     INTEGER iprt
58     real tops(klon),topl(klon),sols(klon),soll(klon)
59     real sens(klon),evap(klon),rain_fall(klon),snow_fall(klon)
60     REAL ts(klon)
61     REAL d_etp_tot, d_qt_tot, d_ec_tot
62     c Output variables
63     REAL fs_bound, fq_bound
64     C
65     C Local variables
66     real stops,stopl,ssols,ssoll
67     real ssens,sfront,slat
68     real airetot, zcpvap, zcwat, zcice
69     REAL rain_fall_tot, snow_fall_tot, evap_tot
70     C
71     integer i
72     C
73     integer pas
74     save pas
75     data pas/0/
76     C
77     pas=pas+1
78     stops=0.
79     stopl=0.
80     ssols=0.
81     ssoll=0.
82     ssens=0.
83     sfront = 0.
84     evap_tot = 0.
85     rain_fall_tot = 0.
86     snow_fall_tot = 0.
87     airetot=0.
88     C
89     C Pour les chaleur specifiques de la vapeur d'eau, de l'eau et de
90     C la glace, on travaille par difference a la chaleur specifique de l'
91     c air sec. En effet, comme on travaille a niveau de pression donne,
92     C toute variation de la masse d'un constituant est totalement
93     c compense par une variation de masse d'air.
94     C
95     zcpvap=RCPV-RCPD
96     zcwat=RCW-RCPD
97     zcice=RCS-RCPD
98     C
99     do i=1,klon
100     stops=stops+tops(i)*airephy(i)
101     stopl=stopl+topl(i)*airephy(i)
102     ssols=ssols+sols(i)*airephy(i)
103     ssoll=ssoll+soll(i)*airephy(i)
104     ssens=ssens+sens(i)*airephy(i)
105     sfront = sfront
106     $ + ( evap(i)*zcpvap-rain_fall(i)*zcwat-snow_fall(i)*zcice
107     $ ) *ts(i) *airephy(i)
108     evap_tot = evap_tot + evap(i)*airephy(i)
109     rain_fall_tot = rain_fall_tot + rain_fall(i)*airephy(i)
110     snow_fall_tot = snow_fall_tot + snow_fall(i)*airephy(i)
111     airetot=airetot+airephy(i)
112     enddo
113     stops=stops/airetot
114     stopl=stopl/airetot
115     ssols=ssols/airetot
116     ssoll=ssoll/airetot
117     ssens=ssens/airetot
118     sfront = sfront/airetot
119     evap_tot = evap_tot /airetot
120     rain_fall_tot = rain_fall_tot/airetot
121     snow_fall_tot = snow_fall_tot/airetot
122     C
123     slat = RLVTT * rain_fall_tot + RLSTT * snow_fall_tot
124     C Heat flux at atm. boundaries
125     fs_bound = stops-stopl - (ssols+ssoll)+ssens+sfront
126     $ + slat
127     C Watter flux at atm. boundaries
128     fq_bound = evap_tot - rain_fall_tot -snow_fall_tot
129     C
130     IF (iprt.ge.1) write(6,6666)
131     $ tit, pas, fs_bound, d_etp_tot, fq_bound, d_qt_tot
132     C
133     IF (iprt.ge.1) write(6,6668)
134     $ tit, pas, d_etp_tot+d_ec_tot-fs_bound, d_qt_tot-fq_bound
135     C
136     IF (iprt.ge.2) write(6,6667)
137     $ tit, pas, stops,stopl,ssols,ssoll,ssens,slat,evap_tot
138     $ ,rain_fall_tot+snow_fall_tot
139    
140     return
141    
142     6666 format('Phys. Flux Budget ',a15,1i6,2f8.2,2(1pE13.5))
143     6667 format('Phys. Boundary Flux ',a15,1i6,6f8.2,2(1pE13.5))
144     6668 format('Phys. Total Budget ',a15,1i6,f8.2,2(1pE13.5))
145    
146     end
147    
148     C======================================================================
149     SUBROUTINE diagetpq(airephy,tit,iprt,idiag,idiag2,dtime
150 guez 10 e ,t,q,ql,qs,u,v,paprs
151 guez 3 s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
152     C======================================================================
153     C
154     C Purpose:
155     C Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
156     C et calcul le flux de chaleur et le flux d'eau necessaire a ces
157     C changements. Ces valeurs sont moyennees sur la surface de tout
158     C le globe et sont exprime en W/2 et kg/s/m2
159     C Outil pour diagnostiquer la conservation de l'energie
160     C et de la masse dans la physique. Suppose que les niveau de
161     c pression entre couche ne varie pas entre 2 appels.
162     C
163     C Plusieurs de ces diagnostics peuvent etre fait en parallele: les
164     c bilans sont sauvegardes dans des tableaux indices. On parlera
165     C "d'indice de diagnostic"
166     c
167     C
168     c======================================================================
169     C Arguments:
170     C airephy-------input-R- grid area
171     C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
172     C iprt----input-I- PRINT level ( <=1 : no PRINT)
173     C idiag---input-I- indice dans lequel sera range les nouveaux
174     C bilans d' entalpie et de masse
175     C idiag2--input-I-les nouveaux bilans d'entalpie et de masse
176     C sont compare au bilan de d'enthalpie de masse de
177     C l'indice numero idiag2
178     C Cas parriculier : si idiag2=0, pas de comparaison, on
179     c sort directement les bilans d'enthalpie et de masse
180     C dtime----input-R- time step (s)
181     c t--------input-R- temperature (K)
182     c q--------input-R- vapeur d'eau (kg/kg)
183     c ql-------input-R- liquid watter (kg/kg)
184     c qs-------input-R- solid watter (kg/kg)
185     c u--------input-R- vitesse u
186     c v--------input-R- vitesse v
187     c paprs----input-R- pression a intercouche (Pa)
188     c
189     C the following total value are computed by UNIT of earth surface
190     C
191     C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
192     c change (J/m2) during one time step (dtime) for the whole
193     C atmosphere (air, watter vapour, liquid and solid)
194     C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
195     C total watter (kg/m2) change during one time step (dtime),
196     C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
197     C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
198     C d_qs------output-R- same, for the solid watter only (kg/m2/s)
199     C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
200     C
201     C other (COMMON...)
202     C RCPD, RCPV, ....
203     C
204     C J.L. Dufresne, July 2002
205     c======================================================================
206    
207     use dimens_m
208     use dimphy
209     use YOMCST
210     use yoethf
211     IMPLICIT NONE
212     C
213     C
214     c Input variables
215     real airephy(klon)
216     CHARACTER*15 tit
217     INTEGER iprt,idiag, idiag2
218     REAL, intent(in):: dtime
219     REAL t(klon,klev), q(klon,klev), ql(klon,klev), qs(klon,klev)
220     REAL u(klon,klev), v(klon,klev)
221     REAL, intent(in):: paprs(klon,klev+1)
222     c Output variables
223     REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
224     C
225     C Local variables
226     c
227     REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
228     . , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
229     c h_vcol_tot-- total enthalpy of vertical air column
230     C (air with watter vapour, liquid and solid) (J/m2)
231     c h_dair_tot-- total enthalpy of dry air (J/m2)
232     c h_qw_tot---- total enthalpy of watter vapour (J/m2)
233     c h_ql_tot---- total enthalpy of liquid watter (J/m2)
234     c h_qs_tot---- total enthalpy of solid watter (J/m2)
235     c qw_tot------ total mass of watter vapour (kg/m2)
236     c ql_tot------ total mass of liquid watter (kg/m2)
237     c qs_tot------ total mass of solid watter (kg/m2)
238     c ec_tot------ total cinetic energy (kg/m2)
239     C
240     REAL zairm(klon,klev) ! layer air mass (kg/m2)
241     REAL zqw_col(klon)
242     REAL zql_col(klon)
243     REAL zqs_col(klon)
244     REAL zec_col(klon)
245     REAL zh_dair_col(klon)
246     REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
247     C
248     REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
249     C
250     REAL airetot, zcpvap, zcwat, zcice
251     C
252     INTEGER i, k
253     C
254     INTEGER ndiag ! max number of diagnostic in parallel
255     PARAMETER (ndiag=10)
256     integer pas(ndiag)
257     save pas
258     data pas/ndiag*0/
259     C
260     REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
261     $ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
262     $ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
263     SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
264     $ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
265    
266     c======================================================================
267     C
268     DO k = 1, klev
269     DO i = 1, klon
270     C layer air mass
271     zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
272     ENDDO
273     END DO
274     C
275     C Reset variables
276     DO i = 1, klon
277     zqw_col(i)=0.
278     zql_col(i)=0.
279     zqs_col(i)=0.
280     zec_col(i) = 0.
281     zh_dair_col(i) = 0.
282     zh_qw_col(i) = 0.
283     zh_ql_col(i) = 0.
284     zh_qs_col(i) = 0.
285     ENDDO
286     C
287     zcpvap=RCPV
288     zcwat=RCW
289     zcice=RCS
290     C
291     C Compute vertical sum for each atmospheric column
292     C ================================================
293     DO k = 1, klev
294     DO i = 1, klon
295     C Watter mass
296     zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)
297     zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)
298     zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)
299     C Cinetic Energy
300     zec_col(i) = zec_col(i)
301     $ +0.5*(u(i,k)**2+v(i,k)**2)*zairm(i,k)
302     C Air enthalpy
303     zh_dair_col(i) = zh_dair_col(i)
304     $ + RCPD*(1.-q(i,k)-ql(i,k)-qs(i,k))*zairm(i,k)*t(i,k)
305     zh_qw_col(i) = zh_qw_col(i)
306     $ + zcpvap*q(i,k)*zairm(i,k)*t(i,k)
307     zh_ql_col(i) = zh_ql_col(i)
308     $ + zcwat*ql(i,k)*zairm(i,k)*t(i,k)
309     $ - RLVTT*ql(i,k)*zairm(i,k)
310     zh_qs_col(i) = zh_qs_col(i)
311     $ + zcice*qs(i,k)*zairm(i,k)*t(i,k)
312     $ - RLSTT*qs(i,k)*zairm(i,k)
313    
314     END DO
315     ENDDO
316     C
317     C Mean over the planete surface
318     C =============================
319     qw_tot = 0.
320     ql_tot = 0.
321     qs_tot = 0.
322     ec_tot = 0.
323     h_vcol_tot = 0.
324     h_dair_tot = 0.
325     h_qw_tot = 0.
326     h_ql_tot = 0.
327     h_qs_tot = 0.
328     airetot=0.
329     C
330     do i=1,klon
331     qw_tot = qw_tot + zqw_col(i)*airephy(i)
332     ql_tot = ql_tot + zql_col(i)*airephy(i)
333     qs_tot = qs_tot + zqs_col(i)*airephy(i)
334     ec_tot = ec_tot + zec_col(i)*airephy(i)
335     h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
336     h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
337     h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
338     h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
339     airetot=airetot+airephy(i)
340     END DO
341     C
342     qw_tot = qw_tot/airetot
343     ql_tot = ql_tot/airetot
344     qs_tot = qs_tot/airetot
345     ec_tot = ec_tot/airetot
346     h_dair_tot = h_dair_tot/airetot
347     h_qw_tot = h_qw_tot/airetot
348     h_ql_tot = h_ql_tot/airetot
349     h_qs_tot = h_qs_tot/airetot
350     C
351     h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
352     C
353     C Compute the change of the atmospheric state compare to the one
354     C stored in "idiag2", and convert it in flux. THis computation
355     C is performed IF idiag2 /= 0 and IF it is not the first CALL
356     c for "idiag"
357     C ===================================
358     C
359     IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
360     d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
361     d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
362     d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
363     d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
364     d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
365     d_qw = (qw_tot - qw_pre(idiag2) )/dtime
366     d_ql = (ql_tot - ql_pre(idiag2) )/dtime
367     d_qs = (qs_tot - qs_pre(idiag2) )/dtime
368     d_ec = (ec_tot - ec_pre(idiag2) )/dtime
369     d_qt = d_qw + d_ql + d_qs
370     ELSE
371     d_h_vcol = 0.
372     d_h_dair = 0.
373     d_h_qw = 0.
374     d_h_ql = 0.
375     d_h_qs = 0.
376     d_qw = 0.
377     d_ql = 0.
378     d_qs = 0.
379     d_ec = 0.
380     d_qt = 0.
381     ENDIF
382     C
383     IF (iprt.ge.2) THEN
384     WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
385     9000 format('Phys. Watter Mass Budget (kg/m2/s)',A15
386     $ ,1i6,10(1pE14.6))
387     WRITE(6,9001) tit,pas(idiag), d_h_vcol
388     9001 format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
389     WRITE(6,9002) tit,pas(idiag), d_ec
390     9002 format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
391     END IF
392     C
393     C Store the new atmospheric state in "idiag"
394     C
395     pas(idiag)=pas(idiag)+1
396     h_vcol_pre(idiag) = h_vcol_tot
397     h_dair_pre(idiag) = h_dair_tot
398     h_qw_pre(idiag) = h_qw_tot
399     h_ql_pre(idiag) = h_ql_tot
400     h_qs_pre(idiag) = h_qs_tot
401     qw_pre(idiag) = qw_tot
402     ql_pre(idiag) = ql_tot
403     qs_pre(idiag) = qs_tot
404     ec_pre (idiag) = ec_tot
405     C
406     RETURN
407     END

  ViewVC Help
Powered by ViewVC 1.1.21