source: LMDZ.3.3/trunk/libf/phylmd/diagphy.F @ 416

Last change on this file since 416 was 374, checked in by lmdzadmin, 22 years ago

Inclusion version initiale
LF

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