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

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

Pllus besoin de comgeomphy
LF

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