source: trunk/LMDZ.TITAN/libf/phytitan/diagphy.F @ 1243

Last change on this file since 1243 was 1048, checked in by slebonnois, 11 years ago

SL: update pour divers details titan + quelques modifs arch et makelmdz

File size: 14.2 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      use cpdet_mod, only: cpdet
212      IMPLICIT NONE
213C
214#include "dimensions.h"
215#include "YOMCST.h"
216C
217c     Input variables
218      real airephy(klon)
219      CHARACTER*15 tit
220      INTEGER iprt,idiag, idiag2
221      REAL dtime
222      REAL t(klon,klev), q(klon,klev), ql(klon,klev), qs(klon,klev)
223      REAL u(klon,klev), v(klon,klev)
224      REAL paprs(klon,klev+1), pplay(klon,klev)
225c     Output variables
226      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
227C
228C     Local variables
229c
230      REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
231     .  , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
232c h_vcol_tot--  total enthalpy of vertical air column
233C            (air with watter vapour, liquid and solid) (J/m2)
234c h_dair_tot-- total enthalpy of dry air (J/m2)
235c h_qw_tot----  total enthalpy of watter vapour (J/m2)
236c h_ql_tot----  total enthalpy of liquid watter (J/m2)
237c h_qs_tot----  total enthalpy of solid watter  (J/m2)
238c qw_tot------  total mass of watter vapour (kg/m2)
239c ql_tot------  total mass of liquid watter (kg/m2)
240c qs_tot------  total mass of solid watter (kg/m2)
241c ec_tot------  total cinetic energy (kg/m2)
242C
243      REAL zairm(klon,klev) ! layer air mass (kg/m2)
244      REAL  zqw_col(klon)
245      REAL  zql_col(klon)
246      REAL  zqs_col(klon)
247      REAL  zec_col(klon)
248      REAL  zh_dair_col(klon)
249      REAL  zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
250C
251      REAL      d_h_dair, d_h_qw, d_h_ql, d_h_qs
252C
253      REAL airetot, zcpvap, zcwat, zcice
254C
255      INTEGER i, k
256C
257      INTEGER ndiag     ! max number of diagnostic in parallel
258      PARAMETER (ndiag=10)
259      integer pas(ndiag)
260      save pas
261      data pas/ndiag*0/
262C     
263      REAL      h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
264     $    , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
265     $    , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
266      SAVE      h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
267     $        , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
268
269c======================================================================
270C
271      DO k = 1, klev
272        DO i = 1, klon
273C         layer air mass
274          zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
275        ENDDO
276      END DO
277C
278C     Reset variables
279      DO i = 1, klon
280        zqw_col(i)=0.
281        zql_col(i)=0.
282        zqs_col(i)=0.
283        zec_col(i) = 0.
284        zh_dair_col(i) = 0.
285        zh_qw_col(i) = 0.
286        zh_ql_col(i) = 0.
287        zh_qs_col(i) = 0.
288      ENDDO
289C
290      zcpvap=RCPV
291      zcwat=RCW
292      zcice=RCS
293C
294C     Compute vertical sum for each atmospheric column
295C     ================================================
296      DO k = 1, klev
297        DO i = 1, klon
298C         Watter mass
299          zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)
300          zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)
301          zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)
302C         Cinetic Energy
303          zec_col(i) =  zec_col(i)
304     $        +0.5*(u(i,k)**2+v(i,k)**2)*zairm(i,k)
305C         Air enthalpy
306! ADAPTATION GCM POUR CP(T)
307          zh_dair_col(i) = zh_dair_col(i)
308     $    + cpdet(t(i,k))*(1.-q(i,k)-ql(i,k)-qs(i,k))*zairm(i,k)*t(i,k)
309          zh_qw_col(i) = zh_qw_col(i)
310     $        + zcpvap*q(i,k)*zairm(i,k)*t(i,k)
311          zh_ql_col(i) = zh_ql_col(i)
312     $        + zcwat*ql(i,k)*zairm(i,k)*t(i,k)
313     $        - RLVTT*ql(i,k)*zairm(i,k)
314          zh_qs_col(i) = zh_qs_col(i)
315     $        + zcice*qs(i,k)*zairm(i,k)*t(i,k)
316     $        - RLSTT*qs(i,k)*zairm(i,k)
317        END DO
318      ENDDO
319C
320C     Mean over the planete surface
321C     =============================
322      qw_tot = 0.
323      ql_tot = 0.
324      qs_tot = 0.
325      ec_tot = 0.
326      h_vcol_tot = 0.
327      h_dair_tot = 0.
328      h_qw_tot = 0.
329      h_ql_tot = 0.
330      h_qs_tot = 0.
331      airetot=0.
332C
333      do i=1,klon
334        qw_tot = qw_tot + zqw_col(i)*airephy(i)
335        ql_tot = ql_tot + zql_col(i)*airephy(i)
336        qs_tot = qs_tot + zqs_col(i)*airephy(i)
337        ec_tot = ec_tot + zec_col(i)*airephy(i)
338        h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
339        h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
340        h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
341        h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
342        airetot=airetot+airephy(i)
343      END DO
344C
345      qw_tot = qw_tot/airetot
346      ql_tot = ql_tot/airetot
347      qs_tot = qs_tot/airetot
348      ec_tot = ec_tot/airetot
349      h_dair_tot = h_dair_tot/airetot
350      h_qw_tot = h_qw_tot/airetot
351      h_ql_tot = h_ql_tot/airetot
352      h_qs_tot = h_qs_tot/airetot
353C
354      h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
355c     print*,'airetot=',airetot,'   h_dair_tot=',h_dair_tot
356C
357C     Compute the change of the atmospheric state compare to the one
358C     stored in "idiag2", and convert it in flux. THis computation
359C     is performed IF idiag2 /= 0 and IF it is not the first CALL
360c     for "idiag"
361C     ===================================
362C
363      IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
364        d_h_vcol  = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
365        d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
366        d_h_qw   = (h_qw_tot  - h_qw_pre(idiag2)  )/dtime
367        d_h_ql   = (h_ql_tot  - h_ql_pre(idiag2)  )/dtime
368        d_h_qs   = (h_qs_tot  - h_qs_pre(idiag2)  )/dtime
369        d_qw     = (qw_tot    - qw_pre(idiag2)    )/dtime
370        d_ql     = (ql_tot    - ql_pre(idiag2)    )/dtime
371        d_qs     = (qs_tot    - qs_pre(idiag2)    )/dtime
372        d_ec     = (ec_tot    - ec_pre(idiag2)    )/dtime
373        d_qt = d_qw + d_ql + d_qs
374      ELSE
375        d_h_vcol = 0.
376        d_h_dair = 0.
377        d_h_qw   = 0.
378        d_h_ql   = 0.
379        d_h_qs   = 0.
380        d_qw     = 0.
381        d_ql     = 0.
382        d_qs     = 0.
383        d_ec     = 0.
384        d_qt     = 0.
385      ENDIF
386C
387      IF (iprt.ge.2) THEN
388        WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
389 9000   format('Phys. Watter Mass Budget (kg/m2/s)',A15
390     $      ,1i6,10(1pE14.6))
391        WRITE(6,9001) tit,pas(idiag), d_h_vcol, h_vcol_tot/dtime
392 9001   format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(E14.6,x))
393        WRITE(6,9002) tit,pas(idiag), d_ec
394 9002   format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F10.2))
395      END IF
396C
397C     Store the new atmospheric state in "idiag"
398C
399      pas(idiag)=pas(idiag)+1
400      h_vcol_pre(idiag)  = h_vcol_tot
401      h_dair_pre(idiag) = h_dair_tot
402      h_qw_pre(idiag)   = h_qw_tot
403      h_ql_pre(idiag)   = h_ql_tot
404      h_qs_pre(idiag)   = h_qs_tot
405      qw_pre(idiag)     = qw_tot
406      ql_pre(idiag)     = ql_tot
407      qs_pre(idiag)     = qs_tot
408      ec_pre (idiag)    = ec_tot
409C
410      RETURN
411      END
Note: See TracBrowser for help on using the repository browser.