source: LMDZ4/trunk/libf/phytherm/diagphy.F @ 1098

Last change on this file since 1098 was 814, checked in by Laurent Fairhead, 17 years ago

Rajout de la physique utilisant les thermiques FH
LF

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