source: lmdz_wrf/branches/LMDZ_WRFmeas/WRFV3/lmdz/diagphy_mod.F90 @ 25

Last change on this file since 25 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.5 KB
Line 
1!
2! $Header$
3!
4MODULE diagphy_mod
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      DO k = 1, klev
281        DO i = 1, klon
282!C         layer air mass
283          zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
284        ENDDO
285      END DO
286!C
287!C     Reset variables
288      DO i = 1, klon
289        zqw_col(i)=0.
290        zql_col(i)=0.
291        zqs_col(i)=0.
292        zec_col(i) = 0.
293        zh_dair_col(i) = 0.
294        zh_qw_col(i) = 0.
295        zh_ql_col(i) = 0.
296        zh_qs_col(i) = 0.
297      ENDDO
298!C
299      zcpvap=RCPV
300      zcwat=RCW
301      zcice=RCS
302!C
303!C     Compute vertical sum for each atmospheric column
304!C     ================================================
305      DO k = 1, klev
306        DO i = 1, klon
307!C         Watter mass
308          zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)
309          zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)
310          zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)
311!C         Cinetic Energy
312          zec_col(i) =  zec_col(i)                                                   &
313       &        +0.5*(u(i,k)**2+v(i,k)**2)*zairm(i,k)
314!C         Air enthalpy
315          zh_dair_col(i) = zh_dair_col(i)                                            &
316       &        + RCPD*(1.-q(i,k)-ql(i,k)-qs(i,k))*zairm(i,k)*t(i,k)
317          zh_qw_col(i) = zh_qw_col(i)                                                &
318       &        + zcpvap*q(i,k)*zairm(i,k)*t(i,k)
319          zh_ql_col(i) = zh_ql_col(i)                                                &
320       &        + zcwat*ql(i,k)*zairm(i,k)*t(i,k)                                    &
321       &        - RLVTT*ql(i,k)*zairm(i,k)
322          zh_qs_col(i) = zh_qs_col(i)                                                &
323       &        + zcice*qs(i,k)*zairm(i,k)*t(i,k)                                    &
324       &        - RLSTT*qs(i,k)*zairm(i,k)
325
326        END DO
327      ENDDO
328
329!C
330!C     Mean over the planete surface
331!C     =============================
332      qw_tot = 0.
333      ql_tot = 0.
334      qs_tot = 0.
335      ec_tot = 0.
336      h_vcol_tot = 0.
337      h_dair_tot = 0.
338      h_qw_tot = 0.
339      h_ql_tot = 0.
340      h_qs_tot = 0.
341      airetot=0.
342!C
343      do i=1,klon
344        qw_tot = qw_tot + zqw_col(i)*airephy(i)
345        ql_tot = ql_tot + zql_col(i)*airephy(i)
346        qs_tot = qs_tot + zqs_col(i)*airephy(i)
347        ec_tot = ec_tot + zec_col(i)*airephy(i)
348        h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
349        h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
350        h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
351        h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
352        airetot=airetot+airephy(i)
353      END DO
354!C
355      qw_tot = qw_tot/airetot
356      ql_tot = ql_tot/airetot
357      qs_tot = qs_tot/airetot
358      ec_tot = ec_tot/airetot
359      h_dair_tot = h_dair_tot/airetot
360      h_qw_tot = h_qw_tot/airetot
361      h_ql_tot = h_ql_tot/airetot
362      h_qs_tot = h_qs_tot/airetot
363!C
364      h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
365!C
366!C     Compute the change of the atmospheric state compare to the one
367!C     stored in "idiag2", and convert it in flux. THis computation
368!C     is performed IF idiag2 /= 0 and IF it is not the first CALL
369!c     for "idiag"
370!C     ===================================
371!C
372      IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
373        d_h_vcol  = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
374        d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
375        d_h_qw   = (h_qw_tot  - h_qw_pre(idiag2)  )/dtime
376        d_h_ql   = (h_ql_tot  - h_ql_pre(idiag2)  )/dtime
377        d_h_qs   = (h_qs_tot  - h_qs_pre(idiag2)  )/dtime
378        d_qw     = (qw_tot    - qw_pre(idiag2)    )/dtime
379        d_ql     = (ql_tot    - ql_pre(idiag2)    )/dtime
380        d_qs     = (qs_tot    - qs_pre(idiag2)    )/dtime
381        d_ec     = (ec_tot    - ec_pre(idiag2)    )/dtime
382        d_qt = d_qw + d_ql + d_qs
383
384      ELSE
385        d_h_vcol = 0.
386        d_h_dair = 0.
387        d_h_qw   = 0.
388        d_h_ql   = 0.
389        d_h_qs   = 0.
390        d_qw     = 0.
391        d_ql     = 0.
392        d_qs     = 0.
393        d_ec     = 0.
394        d_qt     = 0.
395      ENDIF
396!C
397      IF (iprt.ge.2) THEN
398        WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
399 9000   format('Phys. Watter Mass Budget (kg/m2/s)',A15                              &
400       &      ,1i6,10(1pE14.6))
401        WRITE(6,9001) tit,pas(idiag), d_h_vcol
402 9001   format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
403        WRITE(6,9002) tit,pas(idiag), d_ec
404 9002   format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
405      END IF
406!C
407!C     Store the new atmospheric state in "idiag"
408!C
409      pas(idiag)=pas(idiag)+1
410      h_vcol_pre(idiag)  = h_vcol_tot
411      h_dair_pre(idiag) = h_dair_tot
412      h_qw_pre(idiag)   = h_qw_tot
413      h_ql_pre(idiag)   = h_ql_tot
414      h_qs_pre(idiag)   = h_qs_tot
415      qw_pre(idiag)     = qw_tot
416      ql_pre(idiag)     = ql_tot
417      qs_pre(idiag)     = qs_tot
418      ec_pre (idiag)    = ec_tot
419!C
420      RETURN
421      END SUBROUTINE diagetpq
422
423END MODULE diagphy_mod
424
Note: See TracBrowser for help on using the repository browser.