source: LMDZ.3.3/branches/rel-LF/libf/phylmd/diagphy.F @ 436

Last change on this file since 436 was 409, checked in by lmdzadmin, 22 years ago

correction bug JLD
IM/LF

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