source: trunk/libf/phylmd/diagphy.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5

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