source: LMDZ4/trunk/libf/dyn3d/diagedyn.F @ 1381

Last change on this file since 1381 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.7 KB
RevLine 
[524]1!
[1279]2! $Id: diagedyn.F 1279 2009-12-10 09:02:56Z ymeurdesoif $
[524]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"
[1146]60#include "iniprint.h"
[524]61
[1146]62#ifdef CPP_EARTH
[524]63#include "../phylmd/YOMCST.h"
64#include "../phylmd/YOETHF.h"
65#endif
66C
67      INTEGER imjmp1
68      PARAMETER( imjmp1=iim*jjp1)
69c     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
82c     Output variables
83      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
84C
85C     Local variables
86c
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
89c h_vcol_tot--  total enthalpy of vertical air column
90C            (air with watter vapour, liquid and solid) (J/m2)
91c h_dair_tot-- total enthalpy of dry air (J/m2)
92c h_qw_tot----  total enthalpy of watter vapour (J/m2)
93c h_ql_tot----  total enthalpy of liquid watter (J/m2)
94c h_qs_tot----  total enthalpy of solid watter  (J/m2)
95c qw_tot------  total mass of watter vapour (kg/m2)
96c ql_tot------  total mass of liquid watter (kg/m2)
97c qs_tot------  total mass of solid watter (kg/m2)
98c ec_tot------  total cinetic energy (kg/m2)
99C
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)
122C
123      REAL      d_h_dair, d_h_qw, d_h_ql, d_h_qs
124C
125      REAL airetot, zcpvap, zcwat, zcice
126C
127      INTEGER i, k, jj, ij , l ,ip1jjm1
128C
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/
134C     
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
[1146]142#ifdef CPP_EARTH
[524]143c======================================================================
144C     Compute Kinetic enrgy
145      CALL covcont  ( llm    , ucov    , vcov , ucont, vcont        )
146      CALL enercin ( vcov   , ucov  , vcont     , ucont  , ecin  )
147      CALL massdair( p, masse )
148c======================================================================
149C
150C
151      print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?'
152      return
153C     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
162C 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
179C
180C     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
191C
192      zcpvap=RCPV
193      zcwat=RCW
194      zcice=RCS
195C
196C     Compute vertical sum for each atmospheric column
197C     ================================================
198      DO k = 1, llm
199        DO i = 1, imjmp1
200C         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)
204C         Cinetic Energy
205          zec_col(i) =  zec_col(i)
206     $        +zecin(i,k)*zairm(i,k)
207C         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
222C
223C     Mean over the planete surface
224C     =============================
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.
235C
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
247C
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
256C
257      h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
258C
259C     Compute the change of the atmospheric state compare to the one
260C     stored in "idiag2", and convert it in flux. THis computation
261C     is performed IF idiag2 /= 0 and IF it is not the first CALL
262c     for "idiag"
263C     ===================================
264C
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
288C
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))
297C        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
302C
303C     Store the new atmospheric state in "idiag"
304C
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
315C
316#else
[1279]317      write(lunout,*)'diagedyn: Needs Earth physics to function'
[524]318#endif
[1146]319! #endif of #ifdef CPP_EARTH
[524]320      RETURN
321      END
Note: See TracBrowser for help on using the repository browser.