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

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

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

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