source: LMDZ.3.3/trunk/libf/dyn3d/diagedyn.F @ 5306

Last change on this file since 5306 was 415, checked in by lmdzadmin, 22 years ago

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