My Project
 All Classes Files Functions Variables Macros
diagedyn.F
Go to the documentation of this file.
1 !
2 ! $Id: diagedyn.F 1279 2009-12-10 09:02:56Z fairhead $
3 !
4 
5 C======================================================================
6  SUBROUTINE diagedyn(tit,iprt,idiag,idiag2,dtime
7  e , ucov , vcov , ps, p ,pk , teta , q, ql)
8 C======================================================================
9 C
10 C Purpose:
11 C Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
12 C et calcul le flux de chaleur et le flux d'eau necessaire a ces
13 C changements. Ces valeurs sont moyennees sur la surface de tout
14 C le globe et sont exprime en W/2 et kg/s/m2
15 C Outil pour diagnostiquer la conservation de l'energie
16 C et de la masse dans la dynamique.
17 C
18 C
19 c======================================================================
20 C Arguments:
21 C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
22 C iprt----input-I- PRINT level ( <=1 : no PRINT)
23 C idiag---input-I- indice dans lequel sera range les nouveaux
24 C bilans d' entalpie et de masse
25 C idiag2--input-I-les nouveaux bilans d'entalpie et de masse
26 C sont compare au bilan de d'enthalpie de masse de
27 C l'indice numero idiag2
28 C Cas parriculier : si idiag2=0, pas de comparaison, on
29 c sort directement les bilans d'enthalpie et de masse
30 C dtime----input-R- time step (s)
31 C uconv, vconv-input-R- vents covariants (m/s)
32 C ps-------input-R- Surface pressure (Pa)
33 C p--------input-R- pressure at the interfaces
34 C pk-------input-R- pk= (p/Pref)**kappa
35 c teta-----input-R- potential temperature (K)
36 c q--------input-R- vapeur d'eau (kg/kg)
37 c ql-------input-R- liquid watter (kg/kg)
38 c aire-----input-R- mesh surafce (m2)
39 c
40 C the following total value are computed by UNIT of earth surface
41 C
42 C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
43 c change (J/m2) during one time step (dtime) for the whole
44 C atmosphere (air, watter vapour, liquid and solid)
45 C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
46 C total watter (kg/m2) change during one time step (dtime),
47 C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
48 C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
49 C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
50 C
51 C
52 C J.L. Dufresne, July 2002
53 c======================================================================
54 
55  IMPLICIT NONE
56 C
57 #include "dimensions.h"
58 #include "paramet.h"
59 #include "comgeom.h"
60 #include "iniprint.h"
61 
62 #ifdef CPP_EARTH
63 #include "../phylmd/YOMCST.h"
64 #include "../phylmd/YOETHF.h"
65 #endif
66 C
67  INTEGER imjmp1
68  parameter( imjmp1=iim*jjp1)
69 c Input variables
70  CHARACTER*15 tit
71  INTEGER iprt,idiag, idiag2
72  REAL dtime
73  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
74  REAL ps(ip1jmp1) ! pression au sol
75  REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
76  REAL pk (ip1jmp1,llm ) ! = (p/Pref)**kappa
77  REAL teta(ip1jmp1,llm) ! temperature potentielle
78  REAL q(ip1jmp1,llm) ! champs eau vapeur
79  REAL ql(ip1jmp1,llm) ! champs eau liquide
80 
81 
82 c Output variables
83  REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
84 C
85 C Local variables
86 c
87  REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
88  . , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
89 c h_vcol_tot-- total enthalpy of vertical air column
90 C (air with watter vapour, liquid and solid) (J/m2)
91 c h_dair_tot-- total enthalpy of dry air (J/m2)
92 c h_qw_tot---- total enthalpy of watter vapour (J/m2)
93 c h_ql_tot---- total enthalpy of liquid watter (J/m2)
94 c h_qs_tot---- total enthalpy of solid watter (J/m2)
95 c qw_tot------ total mass of watter vapour (kg/m2)
96 c ql_tot------ total mass of liquid watter (kg/m2)
97 c qs_tot------ total mass of solid watter (kg/m2)
98 c ec_tot------ total cinetic energy (kg/m2)
99 C
100  REAL masse(ip1jmp1,llm) ! masse d'air
101  REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
102  REAL ecin(ip1jmp1,llm)
103 
104  REAL zaire(imjmp1)
105  REAL zps(imjmp1)
106  REAL zairm(imjmp1,llm)
107  REAL zecin(imjmp1,llm)
108  REAL zpaprs(imjmp1,llm)
109  REAL zpk(imjmp1,llm)
110  REAL zt(imjmp1,llm)
111  REAL zh(imjmp1,llm)
112  REAL zqw(imjmp1,llm)
113  REAL zql(imjmp1,llm)
114  REAL zqs(imjmp1,llm)
115 
116  REAL zqw_col(imjmp1)
117  REAL zql_col(imjmp1)
118  REAL zqs_col(imjmp1)
119  REAL zec_col(imjmp1)
120  REAL zh_dair_col(imjmp1)
121  REAL zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1)
122 C
123  REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
124 C
125  REAL airetot, zcpvap, zcwat, zcice
126 C
127  INTEGER i, k, jj, ij , l ,ip1jjm1
128 C
129  INTEGER ndiag ! max number of diagnostic in parallel
130  parameter(ndiag=10)
131  integer pas(ndiag)
132  save pas
133  data pas/ndiag*0/
134 C
135  REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
136  $ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
137  $ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
138  SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
139  $ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
140 
141 
142 #ifdef CPP_EARTH
143 c======================================================================
144 C Compute Kinetic enrgy
145  CALL covcont( llm , ucov , vcov , ucont, vcont )
146  CALL enercin( vcov , ucov , vcont , ucont , ecin )
147  CALL massdair( p, masse )
148 c======================================================================
149 C
150 C
151  print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?'
152  return
153 C On ne garde les donnees que dans les colonnes i=1,iim
154  DO jj = 1,jjp1
155  ip1jjm1=iip1*(jj-1)
156  DO ij = 1,iim
157  i=iim*(jj-1)+ij
158  zaire(i)=aire(ij+ip1jjm1)
159  zps(i)=ps(ij+ip1jjm1)
160  ENDDO
161  ENDDO
162 C 3D arrays
163  DO l = 1, llm
164  DO jj = 1,jjp1
165  ip1jjm1=iip1*(jj-1)
166  DO ij = 1,iim
167  i=iim*(jj-1)+ij
168  zairm(i,l) = masse(ij+ip1jjm1,l)
169  zecin(i,l) = ecin(ij+ip1jjm1,l)
170  zpaprs(i,l) = p(ij+ip1jjm1,l)
171  zpk(i,l) = pk(ij+ip1jjm1,l)
172  zh(i,l) = teta(ij+ip1jjm1,l)
173  zqw(i,l) = q(ij+ip1jjm1,l)
174  zql(i,l) = ql(ij+ip1jjm1,l)
175  zqs(i,l) = 0.
176  ENDDO
177  ENDDO
178  ENDDO
179 C
180 C Reset variables
181  DO i = 1, imjmp1
182  zqw_col(i)=0.
183  zql_col(i)=0.
184  zqs_col(i)=0.
185  zec_col(i) = 0.
186  zh_dair_col(i) = 0.
187  zh_qw_col(i) = 0.
188  zh_ql_col(i) = 0.
189  zh_qs_col(i) = 0.
190  ENDDO
191 C
192  zcpvap=rcpv
193  zcwat=rcw
194  zcice=rcs
195 C
196 C Compute vertical sum for each atmospheric column
197 C ================================================
198  DO k = 1, llm
199  DO i = 1, imjmp1
200 C Watter mass
201  zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k)
202  zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k)
203  zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k)
204 C Cinetic Energy
205  zec_col(i) = zec_col(i)
206  $ +zecin(i,k)*zairm(i,k)
207 C Air enthalpy
208  zt(i,k)= zh(i,k) * zpk(i,k) / rcpd
209  zh_dair_col(i) = zh_dair_col(i)
210  $ + rcpd*(1.-zqw(i,k)-zql(i,k)-zqs(i,k))*zairm(i,k)*zt(i,k)
211  zh_qw_col(i) = zh_qw_col(i)
212  $ + zcpvap*zqw(i,k)*zairm(i,k)*zt(i,k)
213  zh_ql_col(i) = zh_ql_col(i)
214  $ + zcwat*zql(i,k)*zairm(i,k)*zt(i,k)
215  $ - rlvtt*zql(i,k)*zairm(i,k)
216  zh_qs_col(i) = zh_qs_col(i)
217  $ + zcice*zqs(i,k)*zairm(i,k)*zt(i,k)
218  $ - rlstt*zqs(i,k)*zairm(i,k)
219 
220  END DO
221  ENDDO
222 C
223 C Mean over the planete surface
224 C =============================
225  qw_tot = 0.
226  ql_tot = 0.
227  qs_tot = 0.
228  ec_tot = 0.
229  h_vcol_tot = 0.
230  h_dair_tot = 0.
231  h_qw_tot = 0.
232  h_ql_tot = 0.
233  h_qs_tot = 0.
234  airetot=0.
235 C
236  do i=1,imjmp1
237  qw_tot = qw_tot + zqw_col(i)
238  ql_tot = ql_tot + zql_col(i)
239  qs_tot = qs_tot + zqs_col(i)
240  ec_tot = ec_tot + zec_col(i)
241  h_dair_tot = h_dair_tot + zh_dair_col(i)
242  h_qw_tot = h_qw_tot + zh_qw_col(i)
243  h_ql_tot = h_ql_tot + zh_ql_col(i)
244  h_qs_tot = h_qs_tot + zh_qs_col(i)
245  airetot=airetot+zaire(i)
246  END DO
247 C
248  qw_tot = qw_tot/airetot
249  ql_tot = ql_tot/airetot
250  qs_tot = qs_tot/airetot
251  ec_tot = ec_tot/airetot
252  h_dair_tot = h_dair_tot/airetot
253  h_qw_tot = h_qw_tot/airetot
254  h_ql_tot = h_ql_tot/airetot
255  h_qs_tot = h_qs_tot/airetot
256 C
257  h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
258 C
259 C Compute the change of the atmospheric state compare to the one
260 C stored in "idiag2", and convert it in flux. THis computation
261 C is performed IF idiag2 /= 0 and IF it is not the first CALL
262 c for "idiag"
263 C ===================================
264 C
265  IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
266  d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
267  d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
268  d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
269  d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
270  d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
271  d_qw = (qw_tot - qw_pre(idiag2) )/dtime
272  d_ql = (ql_tot - ql_pre(idiag2) )/dtime
273  d_qs = (qs_tot - qs_pre(idiag2) )/dtime
274  d_ec = (ec_tot - ec_pre(idiag2) )/dtime
275  d_qt = d_qw + d_ql + d_qs
276  ELSE
277  d_h_vcol = 0.
278  d_h_dair = 0.
279  d_h_qw = 0.
280  d_h_ql = 0.
281  d_h_qs = 0.
282  d_qw = 0.
283  d_ql = 0.
284  d_qs = 0.
285  d_ec = 0.
286  d_qt = 0.
287  ENDIF
288 C
289  IF (iprt.ge.2) THEN
290  WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
291  9000 format('Dyn3d. Watter Mass Budget (kg/m2/s)',a15
292  $ ,1i6,10(1pe14.6))
293  WRITE(6,9001) tit,pas(idiag), d_h_vcol
294  9001 format('Dyn3d. Enthalpy Budget (W/m2) ',a15,1i6,10(f8.2))
295  WRITE(6,9002) tit,pas(idiag), d_ec
296  9002 format('Dyn3d. Cinetic Energy Budget (W/m2) ',a15,1i6,10(f8.2))
297 C WRITE(6,9003) tit,pas(idiag), ec_tot
298  9003 format('Dyn3d. Cinetic Energy (W/m2) ',a15,1i6,10(e15.6))
299  WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec
300  9004 format('Dyn3d. Total Energy Budget (W/m2) ',a15,1i6,10(f8.2))
301  END IF
302 C
303 C Store the new atmospheric state in "idiag"
304 C
305  pas(idiag)=pas(idiag)+1
306  h_vcol_pre(idiag) = h_vcol_tot
307  h_dair_pre(idiag) = h_dair_tot
308  h_qw_pre(idiag) = h_qw_tot
309  h_ql_pre(idiag) = h_ql_tot
310  h_qs_pre(idiag) = h_qs_tot
311  qw_pre(idiag) = qw_tot
312  ql_pre(idiag) = ql_tot
313  qs_pre(idiag) = qs_tot
314  ec_pre(idiag) = ec_tot
315 C
316 #else
317  write(lunout,*)'diagedyn: Needs Earth physics to function'
318 #endif
319 ! #endif of #ifdef CPP_EARTH
320  RETURN
321  END