source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/diagphy.F90 @ 2404

Last change on this file since 2404 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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