source: trunk/libf/phytitan/diagphy.F @ 102

Last change on this file since 102 was 102, checked in by slebonnois, 14 years ago

SL : corrections et modifications dans phytitan correspondant a celles
faites apres compilation Venus. Titan pas encore compile.

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