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

Annotation of /trunk/phylmd/diagphy.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (hide annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 3 months ago) by guez
Original Path: trunk/libf/phylmd/diagphy.f
File size: 14341 byte(s)
Initial import
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     e ,t,q,ql,qs,u,v,paprs,pplay
151     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 pplay----input-R- pression au milieu de couche (Pa)
189     c
190     C the following total value are computed by UNIT of earth surface
191     C
192     C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
193     c change (J/m2) during one time step (dtime) for the whole
194     C atmosphere (air, watter vapour, liquid and solid)
195     C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
196     C total watter (kg/m2) change during one time step (dtime),
197     C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
198     C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
199     C d_qs------output-R- same, for the solid watter only (kg/m2/s)
200     C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
201     C
202     C other (COMMON...)
203     C RCPD, RCPV, ....
204     C
205     C J.L. Dufresne, July 2002
206     c======================================================================
207    
208     use dimens_m
209     use dimphy
210     use YOMCST
211     use yoethf
212     IMPLICIT NONE
213     C
214     C
215     c Input variables
216     real airephy(klon)
217     CHARACTER*15 tit
218     INTEGER iprt,idiag, idiag2
219     REAL, intent(in):: dtime
220     REAL t(klon,klev), q(klon,klev), ql(klon,klev), qs(klon,klev)
221     REAL u(klon,klev), v(klon,klev)
222     REAL, intent(in):: paprs(klon,klev+1)
223     real pplay(klon,klev)
224     c Output variables
225     REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
226     C
227     C Local variables
228     c
229     REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
230     . , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
231     c h_vcol_tot-- total enthalpy of vertical air column
232     C (air with watter vapour, liquid and solid) (J/m2)
233     c h_dair_tot-- total enthalpy of dry air (J/m2)
234     c h_qw_tot---- total enthalpy of watter vapour (J/m2)
235     c h_ql_tot---- total enthalpy of liquid watter (J/m2)
236     c h_qs_tot---- total enthalpy of solid watter (J/m2)
237     c qw_tot------ total mass of watter vapour (kg/m2)
238     c ql_tot------ total mass of liquid watter (kg/m2)
239     c qs_tot------ total mass of solid watter (kg/m2)
240     c ec_tot------ total cinetic energy (kg/m2)
241     C
242     REAL zairm(klon,klev) ! layer air mass (kg/m2)
243     REAL zqw_col(klon)
244     REAL zql_col(klon)
245     REAL zqs_col(klon)
246     REAL zec_col(klon)
247     REAL zh_dair_col(klon)
248     REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
249     C
250     REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
251     C
252     REAL airetot, zcpvap, zcwat, zcice
253     C
254     INTEGER i, k
255     C
256     INTEGER ndiag ! max number of diagnostic in parallel
257     PARAMETER (ndiag=10)
258     integer pas(ndiag)
259     save pas
260     data pas/ndiag*0/
261     C
262     REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
263     $ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
264     $ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
265     SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
266     $ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
267    
268     c======================================================================
269     C
270     DO k = 1, klev
271     DO i = 1, klon
272     C layer air mass
273     zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
274     ENDDO
275     END DO
276     C
277     C Reset variables
278     DO i = 1, klon
279     zqw_col(i)=0.
280     zql_col(i)=0.
281     zqs_col(i)=0.
282     zec_col(i) = 0.
283     zh_dair_col(i) = 0.
284     zh_qw_col(i) = 0.
285     zh_ql_col(i) = 0.
286     zh_qs_col(i) = 0.
287     ENDDO
288     C
289     zcpvap=RCPV
290     zcwat=RCW
291     zcice=RCS
292     C
293     C Compute vertical sum for each atmospheric column
294     C ================================================
295     DO k = 1, klev
296     DO i = 1, klon
297     C Watter mass
298     zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)
299     zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)
300     zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)
301     C Cinetic Energy
302     zec_col(i) = zec_col(i)
303     $ +0.5*(u(i,k)**2+v(i,k)**2)*zairm(i,k)
304     C Air enthalpy
305     zh_dair_col(i) = zh_dair_col(i)
306     $ + RCPD*(1.-q(i,k)-ql(i,k)-qs(i,k))*zairm(i,k)*t(i,k)
307     zh_qw_col(i) = zh_qw_col(i)
308     $ + zcpvap*q(i,k)*zairm(i,k)*t(i,k)
309     zh_ql_col(i) = zh_ql_col(i)
310     $ + zcwat*ql(i,k)*zairm(i,k)*t(i,k)
311     $ - RLVTT*ql(i,k)*zairm(i,k)
312     zh_qs_col(i) = zh_qs_col(i)
313     $ + zcice*qs(i,k)*zairm(i,k)*t(i,k)
314     $ - RLSTT*qs(i,k)*zairm(i,k)
315    
316     END DO
317     ENDDO
318     C
319     C Mean over the planete surface
320     C =============================
321     qw_tot = 0.
322     ql_tot = 0.
323     qs_tot = 0.
324     ec_tot = 0.
325     h_vcol_tot = 0.
326     h_dair_tot = 0.
327     h_qw_tot = 0.
328     h_ql_tot = 0.
329     h_qs_tot = 0.
330     airetot=0.
331     C
332     do i=1,klon
333     qw_tot = qw_tot + zqw_col(i)*airephy(i)
334     ql_tot = ql_tot + zql_col(i)*airephy(i)
335     qs_tot = qs_tot + zqs_col(i)*airephy(i)
336     ec_tot = ec_tot + zec_col(i)*airephy(i)
337     h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
338     h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
339     h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
340     h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
341     airetot=airetot+airephy(i)
342     END DO
343     C
344     qw_tot = qw_tot/airetot
345     ql_tot = ql_tot/airetot
346     qs_tot = qs_tot/airetot
347     ec_tot = ec_tot/airetot
348     h_dair_tot = h_dair_tot/airetot
349     h_qw_tot = h_qw_tot/airetot
350     h_ql_tot = h_ql_tot/airetot
351     h_qs_tot = h_qs_tot/airetot
352     C
353     h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
354     C
355     C Compute the change of the atmospheric state compare to the one
356     C stored in "idiag2", and convert it in flux. THis computation
357     C is performed IF idiag2 /= 0 and IF it is not the first CALL
358     c for "idiag"
359     C ===================================
360     C
361     IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
362     d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
363     d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
364     d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
365     d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
366     d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
367     d_qw = (qw_tot - qw_pre(idiag2) )/dtime
368     d_ql = (ql_tot - ql_pre(idiag2) )/dtime
369     d_qs = (qs_tot - qs_pre(idiag2) )/dtime
370     d_ec = (ec_tot - ec_pre(idiag2) )/dtime
371     d_qt = d_qw + d_ql + d_qs
372     ELSE
373     d_h_vcol = 0.
374     d_h_dair = 0.
375     d_h_qw = 0.
376     d_h_ql = 0.
377     d_h_qs = 0.
378     d_qw = 0.
379     d_ql = 0.
380     d_qs = 0.
381     d_ec = 0.
382     d_qt = 0.
383     ENDIF
384     C
385     IF (iprt.ge.2) THEN
386     WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
387     9000 format('Phys. Watter Mass Budget (kg/m2/s)',A15
388     $ ,1i6,10(1pE14.6))
389     WRITE(6,9001) tit,pas(idiag), d_h_vcol
390     9001 format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
391     WRITE(6,9002) tit,pas(idiag), d_ec
392     9002 format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
393     END IF
394     C
395     C Store the new atmospheric state in "idiag"
396     C
397     pas(idiag)=pas(idiag)+1
398     h_vcol_pre(idiag) = h_vcol_tot
399     h_dair_pre(idiag) = h_dair_tot
400     h_qw_pre(idiag) = h_qw_tot
401     h_ql_pre(idiag) = h_ql_tot
402     h_qs_pre(idiag) = h_qs_tot
403     qw_pre(idiag) = qw_tot
404     ql_pre(idiag) = ql_tot
405     qs_pre(idiag) = qs_tot
406     ec_pre (idiag) = ec_tot
407     C
408     RETURN
409     END

  ViewVC Help
Powered by ViewVC 1.1.21