source: LMDZ4/trunk/libf/dyn3dpar/diagedyn.F @ 1046

Last change on this file since 1046 was 774, checked in by Laurent Fairhead, 17 years ago

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