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

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

Suite de la bascule vers une physique avec thermiques, nouvelle convection, poche froide ...
LF

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