source: lmdz_wrf/trunk/WRFV3/lmdz/diagphy.F90 @ 1465

Last change on this file since 1465 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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