/[lmdze]/trunk/Sources/phylmd/diagetpq.f
ViewVC logotype

Annotation of /trunk/Sources/phylmd/diagetpq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 2 months ago) by guez
Original Path: trunk/phylmd/diagetpq.f
File size: 8509 byte(s)
Changed all ".f90" suffixes to ".f".
1 guez 51 module diagetpq_m
2 guez 47
3 guez 51 IMPLICIT NONE
4 guez 47
5 guez 51 contains
6    
7     SUBROUTINE diagetpq(airephy, tit, iprt, idiag, idiag2, dtime, t, q, ql, qs, &
8     u, v, paprs, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
9    
10     ! From LMDZ4/libf/phylmd/diagphy.F, version 1.1.1.1 2004/05/19 12:53:08
11    
12     ! Purpose:
13    
14     ! Calcule la différence d'enthalpie et de masse d'eau entre deux
15     ! appels et calcule le flux de chaleur et le flux d'eau
16     ! nécessaires à ces changements. Ces valeurs sont moyennées sur la
17     ! surface de tout le globe et sont exprimées en W/m2 et
18     ! kg/s/m2. Outil pour diagnostiquer la conservation de l'énergie
19     ! et de la masse dans la physique. Suppose que les niveaux de
20     ! pression entre les couches ne varient pas entre deux appels.
21    
22     ! Plusieurs de ces diagnostics peuvent être faits en parallèle :
23     ! les bilans sont sauvegardés dans des tableaux indices. On
24     ! parlera "d'indice de diagnostic".
25    
26     ! Jean-Louis Dufresne, July 2002
27    
28 guez 52 USE dimphy, ONLY: klev, klon
29     USE suphec_m, ONLY: rcpd, rcpv, rcs, rcw, rg, rlstt, rlvtt
30 guez 51
31     ! Arguments:
32     ! airephy-------input-R- grid area
33     ! iprt----input-I- PRINT level ( <=1 : no PRINT)
34     ! idiag---input-I- indice dans lequel sera range les nouveaux
35     ! bilans d' entalpie et de masse
36     ! idiag2--input-I-les nouveaux bilans d'entalpie et de masse
37     ! sont compare au bilan de d'enthalpie de masse de
38     ! l'indice numero idiag2
39     ! Cas particulier : si idiag2=0, pas de comparaison, on
40     ! sort directement les bilans d'enthalpie et de masse
41     ! dtime----input-R- time step (s)
42     ! t--------input-R- temperature (K)
43     ! q--------input-R- vapeur d'eau (kg/kg)
44     ! ql-------input-R- liquid water (kg/kg)
45     ! qs-------input-R- solid water (kg/kg)
46     ! u--------input-R- vitesse u
47     ! v--------input-R- vitesse v
48     ! paprs----input-R- pression a intercouche (Pa)
49    
50     ! the following total value are computed by UNIT of earth surface
51    
52     ! d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
53     ! change (J/m2) during one time step (dtime) for the whole
54     ! atmosphere (air, water vapour, liquid and solid)
55     ! d_qt------output-R- total water mass flux (kg/m2/s) defined as the
56     ! total water (kg/m2) change during one time step (dtime),
57     ! d_qw------output-R- same, for the water vapour only (kg/m2/s)
58     ! d_ql------output-R- same, for the liquid water only (kg/m2/s)
59     ! d_qs------output-R- same, for the solid water only (kg/m2/s)
60     ! d_ec------output-R- Kinetic Energy Budget (W/m2) for vertical air column
61    
62     ! other (COMMON...)
63     ! RCPD, RCPV, ....
64    
65     ! Input variables
66     real airephy(klon)
67 guez 62 CHARACTER(len=*), intent(in):: tit ! comment added in PRINT
68 guez 51 INTEGER iprt, idiag, idiag2
69     REAL, intent(in):: dtime
70     REAL, intent(in):: t(klon, klev)
71     REAL, intent(in):: q(klon, klev), ql(klon, klev), qs(klon, klev)
72     REAL u(klon, klev), v(klon, klev)
73     REAL, intent(in):: paprs(klon, klev+1)
74     ! Output variables
75     REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
76    
77     ! Local variables
78    
79     REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot, h_qs_tot, qw_tot, ql_tot
80     real qs_tot , ec_tot
81     ! h_vcol_tot-- total enthalpy of vertical air column
82     ! (air with water vapour, liquid and solid) (J/m2)
83     ! h_dair_tot-- total enthalpy of dry air (J/m2)
84     ! h_qw_tot---- total enthalpy of water vapour (J/m2)
85     ! h_ql_tot---- total enthalpy of liquid water (J/m2)
86     ! h_qs_tot---- total enthalpy of solid water (J/m2)
87     ! qw_tot------ total mass of water vapour (kg/m2)
88     ! ql_tot------ total mass of liquid water (kg/m2)
89     ! qs_tot------ total mass of solid water (kg/m2)
90     ! ec_tot------ total kinetic energy (kg/m2)
91    
92     REAL zairm(klon, klev) ! layer air mass (kg/m2)
93     REAL zqw_col(klon)
94     REAL zql_col(klon)
95     REAL zqs_col(klon)
96     REAL zec_col(klon)
97     REAL zh_dair_col(klon)
98     REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
99    
100     REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
101    
102     REAL airetot, zcpvap, zcwat, zcice
103    
104     INTEGER i, k
105    
106     INTEGER, PARAMETER:: ndiag = 10 ! max number of diagnostic in parallel
107     integer:: pas(ndiag) = 0
108    
109     REAL, save:: h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
110     REAL, save:: h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag), ql_pre(ndiag)
111     REAL, save:: qs_pre(ndiag), ec_pre(ndiag)
112    
113     !-------------------------------------------------------------
114    
115     DO k = 1, klev
116     DO i = 1, klon
117     ! layer air mass
118     zairm(i, k) = (paprs(i, k)-paprs(i, k+1))/RG
119     ENDDO
120     END DO
121    
122     ! Reset variables
123     DO i = 1, klon
124     zqw_col(i)=0.
125     zql_col(i)=0.
126     zqs_col(i)=0.
127     zec_col(i) = 0.
128     zh_dair_col(i) = 0.
129     zh_qw_col(i) = 0.
130     zh_ql_col(i) = 0.
131     zh_qs_col(i) = 0.
132     ENDDO
133    
134     zcpvap=RCPV
135     zcwat=RCW
136     zcice=RCS
137    
138     ! Compute vertical sum for each atmospheric column
139     DO k = 1, klev
140     DO i = 1, klon
141     ! Water mass
142     zqw_col(i) = zqw_col(i) + q(i, k)*zairm(i, k)
143     zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
144     zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
145     ! Kinetic Energy
146     zec_col(i) = zec_col(i) +0.5*(u(i, k)**2+v(i, k)**2)*zairm(i, k)
147     ! Air enthalpy
148     zh_dair_col(i) = zh_dair_col(i) &
149     + RCPD*(1.-q(i, k)-ql(i, k)-qs(i, k))*zairm(i, k)*t(i, k)
150     zh_qw_col(i) = zh_qw_col(i) + zcpvap*q(i, k)*zairm(i, k)*t(i, k)
151     zh_ql_col(i) = zh_ql_col(i) &
152     + zcwat*ql(i, k)*zairm(i, k)*t(i, k) &
153     - RLVTT*ql(i, k)*zairm(i, k)
154     zh_qs_col(i) = zh_qs_col(i) &
155     + zcice*qs(i, k)*zairm(i, k)*t(i, k) &
156     - RLSTT*qs(i, k)*zairm(i, k)
157     END DO
158     ENDDO
159    
160     ! Mean over the planet surface
161     qw_tot = 0.
162     ql_tot = 0.
163     qs_tot = 0.
164     ec_tot = 0.
165     h_vcol_tot = 0.
166     h_dair_tot = 0.
167     h_qw_tot = 0.
168     h_ql_tot = 0.
169     h_qs_tot = 0.
170     airetot=0.
171    
172     do i=1, klon
173     qw_tot = qw_tot + zqw_col(i)*airephy(i)
174     ql_tot = ql_tot + zql_col(i)*airephy(i)
175     qs_tot = qs_tot + zqs_col(i)*airephy(i)
176     ec_tot = ec_tot + zec_col(i)*airephy(i)
177     h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
178     h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
179     h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
180     h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
181     airetot=airetot+airephy(i)
182     END DO
183    
184     qw_tot = qw_tot/airetot
185     ql_tot = ql_tot/airetot
186     qs_tot = qs_tot/airetot
187     ec_tot = ec_tot/airetot
188     h_dair_tot = h_dair_tot/airetot
189     h_qw_tot = h_qw_tot/airetot
190     h_ql_tot = h_ql_tot/airetot
191     h_qs_tot = h_qs_tot/airetot
192    
193     h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
194    
195     ! Compute the change of the atmospheric state compared to the one
196     ! stored in "idiag2", and convert it in flux. This computation is
197     ! performed if idiag2 /= 0 and if it is not the first call for
198     ! "idiag".
199    
200     IF ((idiag2 > 0) .and. (pas(idiag2) /= 0)) THEN
201     d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
202     d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
203     d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
204     d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
205     d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
206     d_qw = (qw_tot - qw_pre(idiag2) )/dtime
207     d_ql = (ql_tot - ql_pre(idiag2) )/dtime
208     d_qs = (qs_tot - qs_pre(idiag2) )/dtime
209     d_ec = (ec_tot - ec_pre(idiag2) )/dtime
210     d_qt = d_qw + d_ql + d_qs
211     ELSE
212     d_h_vcol = 0.
213     d_h_dair = 0.
214     d_h_qw = 0.
215     d_h_ql = 0.
216     d_h_qs = 0.
217     d_qw = 0.
218     d_ql = 0.
219     d_qs = 0.
220     d_ec = 0.
221     d_qt = 0.
222     ENDIF
223    
224     IF (iprt >= 2) THEN
225     WRITE(6, 9000) tit, pas(idiag), d_qt, d_qw, d_ql, d_qs
226     9000 format('Phys. Water Mass Budget (kg/m2/s)', A15, 1i6, 10(1pE14.6))
227     WRITE(6, 9001) tit, pas(idiag), d_h_vcol
228     9001 format('Phys. Enthalpy Budget (W/m2) ', A15, 1i6, 10(F8.2))
229     WRITE(6, 9002) tit, pas(idiag), d_ec
230     9002 format('Phys. Kinetic Energy Budget (W/m2) ', A15, 1i6, 10(F8.2))
231     END IF
232    
233     ! Store the new atmospheric state in "idiag"
234     pas(idiag)=pas(idiag)+1
235     h_vcol_pre(idiag) = h_vcol_tot
236     h_dair_pre(idiag) = h_dair_tot
237     h_qw_pre(idiag) = h_qw_tot
238     h_ql_pre(idiag) = h_ql_tot
239     h_qs_pre(idiag) = h_qs_tot
240     qw_pre(idiag) = qw_tot
241     ql_pre(idiag) = ql_tot
242     qs_pre(idiag) = qs_tot
243     ec_pre (idiag) = ec_tot
244    
245     END SUBROUTINE diagetpq
246    
247     end module diagetpq_m

  ViewVC Help
Powered by ViewVC 1.1.21