source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/diagedyn.F @ 5444

Last change on this file since 5444 was 545, checked in by lmdzadmin, 20 years ago

Suite des modifs pour adaptation a Prism MED
LF
  • 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"
[545]57Cmed#include "../phylmd/YOMCST.h"
58Cmed#include "../phylmd/YOETHF.h"
59#include "YOMCST.h"
60#include "YOETHF.h"
[415]61C
62      INTEGER imjmp1
63      PARAMETER( imjmp1=iim*jjp1)
64c     Input variables
65      CHARACTER*15 tit
66      INTEGER iprt,idiag, idiag2
67      REAL dtime
68      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
69      REAL ps(ip1jmp1)                       ! pression  au sol
70      REAL p (ip1jmp1,llmp1  )  ! pression aux interfac.des couches
[431]71      REAL pk (ip1jmp1,llm  )  ! = (p/Pref)**kappa
[415]72      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
73      REAL q(ip1jmp1,llm)               ! champs eau vapeur
74      REAL ql(ip1jmp1,llm)               ! champs eau liquide
75      REAL aire(ip1jmp1)               ! aire des mailles
76
77
78c     Output variables
79      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
80C
81C     Local variables
82c
83      REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
84     .  , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
85c h_vcol_tot--  total enthalpy of vertical air column
86C            (air with watter vapour, liquid and solid) (J/m2)
87c h_dair_tot-- total enthalpy of dry air (J/m2)
88c h_qw_tot----  total enthalpy of watter vapour (J/m2)
89c h_ql_tot----  total enthalpy of liquid watter (J/m2)
90c h_qs_tot----  total enthalpy of solid watter  (J/m2)
91c qw_tot------  total mass of watter vapour (kg/m2)
92c ql_tot------  total mass of liquid watter (kg/m2)
93c qs_tot------  total mass of solid watter (kg/m2)
94c ec_tot------  total cinetic energy (kg/m2)
95C
96      REAL masse(ip1jmp1,llm)                ! masse d'air
97      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
98      REAL ecin(ip1jmp1,llm)
99
100      REAL zaire(imjmp1)
101      REAL zps(imjmp1)
102      REAL zairm(imjmp1,llm)
103      REAL zecin(imjmp1,llm)
104      REAL zpaprs(imjmp1,llm)
105      REAL zpk(imjmp1,llm)
106      REAL zt(imjmp1,llm)
107      REAL zh(imjmp1,llm)
108      REAL zqw(imjmp1,llm)
109      REAL zql(imjmp1,llm)
110      REAL zqs(imjmp1,llm)
111
112      REAL  zqw_col(imjmp1)
113      REAL  zql_col(imjmp1)
114      REAL  zqs_col(imjmp1)
115      REAL  zec_col(imjmp1)
116      REAL  zh_dair_col(imjmp1)
117      REAL  zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1)
118C
119      REAL      d_h_dair, d_h_qw, d_h_ql, d_h_qs
120C
121      REAL airetot, zcpvap, zcwat, zcice
122C
123      INTEGER i, k, jj, ij , l ,ip1jjm1
124C
125      INTEGER ndiag     ! max number of diagnostic in parallel
126      PARAMETER (ndiag=10)
127      integer pas(ndiag)
128      save pas
129      data pas/ndiag*0/
130C     
131      REAL      h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
132     $    , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
133     $    , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
134      SAVE      h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
135     $        , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
136
137
138c======================================================================
139C     Compute Kinetic enrgy
140      CALL covcont  ( llm    , ucov    , vcov , ucont, vcont        )
141      CALL enercin ( vcov   , ucov  , vcont     , ucont  , ecin  )
142      CALL massdair( p, masse )
143c======================================================================
144C
145C
146C     On ne garde les donnees que dans les colonnes i=1,iim
147      DO jj = 1,jjp1
148        ip1jjm1=iip1*(jj-1)
149        DO ij =  1,iim
150          i=iim*(jj-1)+ij
151          zaire(i)=aire(ij+ip1jjm1)
152          zps(i)=ps(ij+ip1jjm1)
153        ENDDO
154      ENDDO
155C 3D arrays
156      DO l  =  1, llm
157        DO jj = 1,jjp1
158          ip1jjm1=iip1*(jj-1)
159          DO ij =  1,iim
160            i=iim*(jj-1)+ij
161            zairm(i,l) = masse(ij+ip1jjm1,l)
162            zecin(i,l) = ecin(ij+ip1jjm1,l)
163            zpaprs(i,l) = p(ij+ip1jjm1,l)
164            zpk(i,l) = pk(ij+ip1jjm1,l)
165            zh(i,l) = teta(ij+ip1jjm1,l)
166            zqw(i,l) = q(ij+ip1jjm1,l)
167            zql(i,l) = ql(ij+ip1jjm1,l)
168            zqs(i,l) = 0.
169          ENDDO
170        ENDDO
171      ENDDO
172C
173C     Reset variables
174      DO i = 1, imjmp1
175        zqw_col(i)=0.
176        zql_col(i)=0.
177        zqs_col(i)=0.
178        zec_col(i) = 0.
179        zh_dair_col(i) = 0.
180        zh_qw_col(i) = 0.
181        zh_ql_col(i) = 0.
182        zh_qs_col(i) = 0.
183      ENDDO
184C
185      zcpvap=RCPV
186      zcwat=RCW
187      zcice=RCS
188C
189C     Compute vertical sum for each atmospheric column
190C     ================================================
191      DO k = 1, llm
192        DO i = 1, imjmp1
193C         Watter mass
194          zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k)
195          zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k)
196          zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k)
197C         Cinetic Energy
198          zec_col(i) =  zec_col(i)
199     $        +zecin(i,k)*zairm(i,k)
200C         Air enthalpy
201          zt(i,k)= zh(i,k) * zpk(i,k) / RCPD
202          zh_dair_col(i) = zh_dair_col(i)
203     $        + RCPD*(1.-zqw(i,k)-zql(i,k)-zqs(i,k))*zairm(i,k)*zt(i,k)
204          zh_qw_col(i) = zh_qw_col(i)
205     $        + zcpvap*zqw(i,k)*zairm(i,k)*zt(i,k)
206          zh_ql_col(i) = zh_ql_col(i)
207     $        + zcwat*zql(i,k)*zairm(i,k)*zt(i,k)
208     $        - RLVTT*zql(i,k)*zairm(i,k)
209          zh_qs_col(i) = zh_qs_col(i)
210     $        + zcice*zqs(i,k)*zairm(i,k)*zt(i,k)
211     $        - RLSTT*zqs(i,k)*zairm(i,k)
212
213        END DO
214      ENDDO
215C
216C     Mean over the planete surface
217C     =============================
218      qw_tot = 0.
219      ql_tot = 0.
220      qs_tot = 0.
221      ec_tot = 0.
222      h_vcol_tot = 0.
223      h_dair_tot = 0.
224      h_qw_tot = 0.
225      h_ql_tot = 0.
226      h_qs_tot = 0.
227      airetot=0.
228C
229      do i=1,imjmp1
230        qw_tot = qw_tot + zqw_col(i)
231        ql_tot = ql_tot + zql_col(i)
232        qs_tot = qs_tot + zqs_col(i)
233        ec_tot = ec_tot + zec_col(i)
234        h_dair_tot = h_dair_tot + zh_dair_col(i)
235        h_qw_tot = h_qw_tot + zh_qw_col(i)
236        h_ql_tot = h_ql_tot + zh_ql_col(i)
237        h_qs_tot = h_qs_tot + zh_qs_col(i)
238        airetot=airetot+zaire(i)
239      END DO
240C
241      qw_tot = qw_tot/airetot
242      ql_tot = ql_tot/airetot
243      qs_tot = qs_tot/airetot
244      ec_tot = ec_tot/airetot
245      h_dair_tot = h_dair_tot/airetot
246      h_qw_tot = h_qw_tot/airetot
247      h_ql_tot = h_ql_tot/airetot
248      h_qs_tot = h_qs_tot/airetot
249C
250      h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
251C
252C     Compute the change of the atmospheric state compare to the one
253C     stored in "idiag2", and convert it in flux. THis computation
254C     is performed IF idiag2 /= 0 and IF it is not the first CALL
255c     for "idiag"
256C     ===================================
257C
258      IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
259        d_h_vcol  = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
260        d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
261        d_h_qw   = (h_qw_tot  - h_qw_pre(idiag2)  )/dtime
262        d_h_ql   = (h_ql_tot  - h_ql_pre(idiag2)  )/dtime
263        d_h_qs   = (h_qs_tot  - h_qs_pre(idiag2)  )/dtime
264        d_qw     = (qw_tot    - qw_pre(idiag2)    )/dtime
265        d_ql     = (ql_tot    - ql_pre(idiag2)    )/dtime
266        d_qs     = (qs_tot    - qs_pre(idiag2)    )/dtime
267        d_ec     = (ec_tot    - ec_pre(idiag2)    )/dtime
268        d_qt = d_qw + d_ql + d_qs
269      ELSE
270        d_h_vcol = 0.
271        d_h_dair = 0.
272        d_h_qw   = 0.
273        d_h_ql   = 0.
274        d_h_qs   = 0.
275        d_qw     = 0.
276        d_ql     = 0.
277        d_qs     = 0.
278        d_ec     = 0.
279        d_qt     = 0.
280      ENDIF
281C
282      IF (iprt.ge.2) THEN
283        WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
284 9000   format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15
285     $      ,1i6,10(1pE14.6))
286        WRITE(6,9001) tit,pas(idiag), d_h_vcol
287 9001   format('Dyn3d. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
288        WRITE(6,9002) tit,pas(idiag), d_ec
289 9002   format('Dyn3d. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
290C        WRITE(6,9003) tit,pas(idiag), ec_tot
291 9003   format('Dyn3d. Cinetic Energy (W/m2) ',A15,1i6,10(E15.6))
292        WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec
293 9004   format('Dyn3d. Total Energy Budget (W/m2) ',A15,1i6,10(F8.2))
294      END IF
295C
296C     Store the new atmospheric state in "idiag"
297C
298      pas(idiag)=pas(idiag)+1
299      h_vcol_pre(idiag)  = h_vcol_tot
300      h_dair_pre(idiag) = h_dair_tot
301      h_qw_pre(idiag)   = h_qw_tot
302      h_ql_pre(idiag)   = h_ql_tot
303      h_qs_pre(idiag)   = h_qs_tot
304      qw_pre(idiag)     = qw_tot
305      ql_pre(idiag)     = ql_tot
306      qs_pre(idiag)     = qs_tot
307      ec_pre (idiag)    = ec_tot
308C
309      RETURN
310      END
Note: See TracBrowser for help on using the repository browser.