source: lmdz_wrf/trunk/WRFV3/lmdz/diagphy_mod.F90 @ 31

Last change on this file since 31 was 31, checked in by lfita, 10 years ago

Adding checking for NaN

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