source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/diagedyn.F @ 5097

Last change on this file since 5097 was 5089, checked in by abarral, 4 months ago

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