source: LMDZ4/trunk/libf/phylmd/diagphy.F @ 602

Last change on this file since 602 was 524, checked in by lmdzadmin, 20 years ago

Initial revision

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