source: LMDZ6/trunk/libf/dyn3d_common/diagedyn.F @ 4651

Last change on this file since 4651 was 4593, checked in by yann meurdesoif, 18 months ago

Replace #include (c preprocessor) by INCLUDE (fortran keyword)

in phylmd (except rrtm and ecrad) filtrez, dy3dmem and dyn3dcommon

Other directories will follow
YM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.3 KB
RevLine 
[524]1!
[1279]2! $Id: diagedyn.F 4593 2023-06-29 13:55:54Z fhourdin $
[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 
[2239]55      USE control_mod, ONLY : planet_type
56     
[524]57      IMPLICIT NONE
58C
[4593]59      INCLUDE "dimensions.h"
60      INCLUDE "paramet.h"
61      INCLUDE "comgeom.h"
62      INCLUDE "iniprint.h"
[524]63
[2239]64!#ifdef CPP_EARTH
[4593]65!      INCLUDE "../phylmd/YOMCST.h"
66!      INCLUDE "../phylmd/YOETHF.h"
[2239]67!#endif
68! Ehouarn: for now set these parameters to what is in Earth physics...
69!          (cf ../phylmd/suphel.h)
70!          this should be generalized...
71      REAL,PARAMETER :: RCPD=
72     &               3.5*(1000.*(6.0221367E+23*1.380658E-23)/28.9644)
73      REAL,PARAMETER :: RCPV=
74     &               4.*(1000.*(6.0221367E+23*1.380658E-23)/18.0153)
75      REAL,PARAMETER :: RCS=RCPV
76      REAL,PARAMETER :: RCW=RCPV
77      REAL,PARAMETER :: RLSTT=2.8345E+6
78      REAL,PARAMETER :: RLVTT=2.5008E+6
79!
[524]80C
81      INTEGER imjmp1
82      PARAMETER( imjmp1=iim*jjp1)
83c     Input variables
84      CHARACTER*15 tit
85      INTEGER iprt,idiag, idiag2
86      REAL dtime
87      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
88      REAL ps(ip1jmp1)                       ! pression  au sol
89      REAL p (ip1jmp1,llmp1  )  ! pression aux interfac.des couches
90      REAL pk (ip1jmp1,llm  )  ! = (p/Pref)**kappa
91      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
92      REAL q(ip1jmp1,llm)               ! champs eau vapeur
93      REAL ql(ip1jmp1,llm)               ! champs eau liquide
94
95
96c     Output variables
97      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
98C
99C     Local variables
100c
101      REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
102     .  , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
103c h_vcol_tot--  total enthalpy of vertical air column
104C            (air with watter vapour, liquid and solid) (J/m2)
105c h_dair_tot-- total enthalpy of dry air (J/m2)
106c h_qw_tot----  total enthalpy of watter vapour (J/m2)
107c h_ql_tot----  total enthalpy of liquid watter (J/m2)
108c h_qs_tot----  total enthalpy of solid watter  (J/m2)
109c qw_tot------  total mass of watter vapour (kg/m2)
110c ql_tot------  total mass of liquid watter (kg/m2)
111c qs_tot------  total mass of solid watter (kg/m2)
112c ec_tot------  total cinetic energy (kg/m2)
113C
114      REAL masse(ip1jmp1,llm)                ! masse d'air
115      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
116      REAL ecin(ip1jmp1,llm)
117
118      REAL zaire(imjmp1)
119      REAL zps(imjmp1)
120      REAL zairm(imjmp1,llm)
121      REAL zecin(imjmp1,llm)
122      REAL zpaprs(imjmp1,llm)
123      REAL zpk(imjmp1,llm)
124      REAL zt(imjmp1,llm)
125      REAL zh(imjmp1,llm)
126      REAL zqw(imjmp1,llm)
127      REAL zql(imjmp1,llm)
128      REAL zqs(imjmp1,llm)
129
130      REAL  zqw_col(imjmp1)
131      REAL  zql_col(imjmp1)
132      REAL  zqs_col(imjmp1)
133      REAL  zec_col(imjmp1)
134      REAL  zh_dair_col(imjmp1)
135      REAL  zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1)
136C
137      REAL      d_h_dair, d_h_qw, d_h_ql, d_h_qs
138C
139      REAL airetot, zcpvap, zcwat, zcice
140C
141      INTEGER i, k, jj, ij , l ,ip1jjm1
142C
143      INTEGER ndiag     ! max number of diagnostic in parallel
144      PARAMETER (ndiag=10)
145      integer pas(ndiag)
146      save pas
147      data pas/ndiag*0/
148C     
149      REAL      h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
150     $    , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
151     $    , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
152      SAVE      h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
153     $        , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
154
155
[2239]156!#ifdef CPP_EARTH
157      IF (planet_type=="earth") THEN
158     
[524]159c======================================================================
160C     Compute Kinetic enrgy
161      CALL covcont  ( llm    , ucov    , vcov , ucont, vcont        )
162      CALL enercin ( vcov   , ucov  , vcont     , ucont  , ecin  )
163      CALL massdair( p, masse )
164c======================================================================
165C
166C
167      print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?'
168      return
169C     On ne garde les donnees que dans les colonnes i=1,iim
170      DO jj = 1,jjp1
171        ip1jjm1=iip1*(jj-1)
172        DO ij =  1,iim
173          i=iim*(jj-1)+ij
174          zaire(i)=aire(ij+ip1jjm1)
175          zps(i)=ps(ij+ip1jjm1)
176        ENDDO
177      ENDDO
178C 3D arrays
179      DO l  =  1, llm
180        DO jj = 1,jjp1
181          ip1jjm1=iip1*(jj-1)
182          DO ij =  1,iim
183            i=iim*(jj-1)+ij
184            zairm(i,l) = masse(ij+ip1jjm1,l)
185            zecin(i,l) = ecin(ij+ip1jjm1,l)
186            zpaprs(i,l) = p(ij+ip1jjm1,l)
187            zpk(i,l) = pk(ij+ip1jjm1,l)
188            zh(i,l) = teta(ij+ip1jjm1,l)
189            zqw(i,l) = q(ij+ip1jjm1,l)
190            zql(i,l) = ql(ij+ip1jjm1,l)
191            zqs(i,l) = 0.
192          ENDDO
193        ENDDO
194      ENDDO
195C
196C     Reset variables
197      DO i = 1, imjmp1
198        zqw_col(i)=0.
199        zql_col(i)=0.
200        zqs_col(i)=0.
201        zec_col(i) = 0.
202        zh_dair_col(i) = 0.
203        zh_qw_col(i) = 0.
204        zh_ql_col(i) = 0.
205        zh_qs_col(i) = 0.
206      ENDDO
207C
208      zcpvap=RCPV
209      zcwat=RCW
210      zcice=RCS
211C
212C     Compute vertical sum for each atmospheric column
213C     ================================================
214      DO k = 1, llm
215        DO i = 1, imjmp1
216C         Watter mass
217          zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k)
218          zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k)
219          zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k)
220C         Cinetic Energy
221          zec_col(i) =  zec_col(i)
222     $        +zecin(i,k)*zairm(i,k)
223C         Air enthalpy
224          zt(i,k)= zh(i,k) * zpk(i,k) / RCPD
225          zh_dair_col(i) = zh_dair_col(i)
226     $        + RCPD*(1.-zqw(i,k)-zql(i,k)-zqs(i,k))*zairm(i,k)*zt(i,k)
227          zh_qw_col(i) = zh_qw_col(i)
228     $        + zcpvap*zqw(i,k)*zairm(i,k)*zt(i,k)
229          zh_ql_col(i) = zh_ql_col(i)
230     $        + zcwat*zql(i,k)*zairm(i,k)*zt(i,k)
231     $        - RLVTT*zql(i,k)*zairm(i,k)
232          zh_qs_col(i) = zh_qs_col(i)
233     $        + zcice*zqs(i,k)*zairm(i,k)*zt(i,k)
234     $        - RLSTT*zqs(i,k)*zairm(i,k)
235
236        END DO
237      ENDDO
238C
239C     Mean over the planete surface
240C     =============================
241      qw_tot = 0.
242      ql_tot = 0.
243      qs_tot = 0.
244      ec_tot = 0.
245      h_vcol_tot = 0.
246      h_dair_tot = 0.
247      h_qw_tot = 0.
248      h_ql_tot = 0.
249      h_qs_tot = 0.
250      airetot=0.
251C
252      do i=1,imjmp1
253        qw_tot = qw_tot + zqw_col(i)
254        ql_tot = ql_tot + zql_col(i)
255        qs_tot = qs_tot + zqs_col(i)
256        ec_tot = ec_tot + zec_col(i)
257        h_dair_tot = h_dair_tot + zh_dair_col(i)
258        h_qw_tot = h_qw_tot + zh_qw_col(i)
259        h_ql_tot = h_ql_tot + zh_ql_col(i)
260        h_qs_tot = h_qs_tot + zh_qs_col(i)
261        airetot=airetot+zaire(i)
262      END DO
263C
264      qw_tot = qw_tot/airetot
265      ql_tot = ql_tot/airetot
266      qs_tot = qs_tot/airetot
267      ec_tot = ec_tot/airetot
268      h_dair_tot = h_dair_tot/airetot
269      h_qw_tot = h_qw_tot/airetot
270      h_ql_tot = h_ql_tot/airetot
271      h_qs_tot = h_qs_tot/airetot
272C
273      h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
274C
275C     Compute the change of the atmospheric state compare to the one
276C     stored in "idiag2", and convert it in flux. THis computation
277C     is performed IF idiag2 /= 0 and IF it is not the first CALL
278c     for "idiag"
279C     ===================================
280C
281      IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
282        d_h_vcol  = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
283        d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
284        d_h_qw   = (h_qw_tot  - h_qw_pre(idiag2)  )/dtime
285        d_h_ql   = (h_ql_tot  - h_ql_pre(idiag2)  )/dtime
286        d_h_qs   = (h_qs_tot  - h_qs_pre(idiag2)  )/dtime
287        d_qw     = (qw_tot    - qw_pre(idiag2)    )/dtime
288        d_ql     = (ql_tot    - ql_pre(idiag2)    )/dtime
289        d_qs     = (qs_tot    - qs_pre(idiag2)    )/dtime
290        d_ec     = (ec_tot    - ec_pre(idiag2)    )/dtime
291        d_qt = d_qw + d_ql + d_qs
292      ELSE
293        d_h_vcol = 0.
294        d_h_dair = 0.
295        d_h_qw   = 0.
296        d_h_ql   = 0.
297        d_h_qs   = 0.
298        d_qw     = 0.
299        d_ql     = 0.
300        d_qs     = 0.
301        d_ec     = 0.
302        d_qt     = 0.
303      ENDIF
304C
305      IF (iprt.ge.2) THEN
306        WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
307 9000   format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15
308     $      ,1i6,10(1pE14.6))
309        WRITE(6,9001) tit,pas(idiag), d_h_vcol
310 9001   format('Dyn3d. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
311        WRITE(6,9002) tit,pas(idiag), d_ec
312 9002   format('Dyn3d. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
313C        WRITE(6,9003) tit,pas(idiag), ec_tot
314 9003   format('Dyn3d. Cinetic Energy (W/m2) ',A15,1i6,10(E15.6))
315        WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec
316 9004   format('Dyn3d. Total Energy Budget (W/m2) ',A15,1i6,10(F8.2))
317      END IF
318C
319C     Store the new atmospheric state in "idiag"
320C
321      pas(idiag)=pas(idiag)+1
322      h_vcol_pre(idiag)  = h_vcol_tot
323      h_dair_pre(idiag) = h_dair_tot
324      h_qw_pre(idiag)   = h_qw_tot
325      h_ql_pre(idiag)   = h_ql_tot
326      h_qs_pre(idiag)   = h_qs_tot
327      qw_pre(idiag)     = qw_tot
328      ql_pre(idiag)     = ql_tot
329      qs_pre(idiag)     = qs_tot
330      ec_pre (idiag)    = ec_tot
331C
[2239]332!#else
333      ELSE
334        write(lunout,*)'diagedyn: set to function with Earth parameters'
335      ENDIF ! of if (planet_type=="earth")
336!#endif
[1146]337! #endif of #ifdef CPP_EARTH
[524]338      RETURN
339      END
Note: See TracBrowser for help on using the repository browser.