source: trunk/LMDZ.VENUS/libf/phyvenus/diagphy.F @ 1723

Last change on this file since 1723 was 1621, checked in by emillour, 8 years ago

Further work on full dynamics/physics separation.

LMDZ.COMMON:

  • added phy_common/vertical_layers_mod.F90 to store information on vertical grid. This is where routines in the physics should get the information.
  • The contents of vertical_layers_mod intialized via dynphy_lonlat/inigeomphy_mod.F90.

LMDZ.MARS:

  • physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
  • created an "ini_tracer_mod" routine in module "tracer_mod" for a cleaner initialization of the later.
  • removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.GENERIC:

  • physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
  • added nqtot to tracer_h.F90.
  • removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.VENUS:

  • physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
  • added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
  • added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics. Initialized via iniphysiq. IMPORTANT: there are some hard-coded constants! These should match what is in cpdet_mod.F90 in the dynamics.
  • got rid of references to moyzon_mod module within the physics. The required variables (tmoy, plevmoy) are passed to the physics as arguments to physiq.

LMDZ.TITAN:

  • added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
  • added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics.
  • Extra work required to completely decouple physics and dynamics: moyzon_mod should be cleaned up and information passed from dynamics to physics as as arguments. Likewise moyzon_ch and moyzon_mu should not be queried from logic_mod (which is in the dynamics).

EM
File size: 14.1 KB
Line 
1!
2! $Header: /home/cvsroot/LMDZ4/libf/phylmd/diagphy.F,v 1.1.1.1 2004/05/19 12:53:08 lmdzadmin Exp $
3!
4      SUBROUTINE diagphy(airephy,tit,iprt
5     $    , tops, topl, sols, soll, sens
6     $    , evap, rain_fall, snow_fall, ts
7     $    , d_etp_tot, d_qt_tot, d_ec_tot
8     $    , fs_bound, fq_bound)
9
10! ATTENTION !! PAS DU TOUT A JOUR POUR VENUS OU TITAN...
11
12C======================================================================
13C
14C Purpose:
15C    Compute the thermal flux and the watter mass flux at the atmosphere
16c    boundaries. Print them and also the atmospheric enthalpy change and
17C    the  atmospheric mass change.
18C
19C Arguments:
20C airephy-------input-R-  grid area
21C tit---------input-A15- Comment to be added in PRINT (CHARACTER*15)
22C iprt--------input-I-  PRINT level ( <=0 : no PRINT)
23C tops(klon)--input-R-  SW rad. at TOA (W/m2), positive up.
24C topl(klon)--input-R-  LW rad. at TOA (W/m2), positive down
25C sols(klon)--input-R-  Net SW flux above surface (W/m2), positive up
26C                   (i.e. -1 * flux absorbed by the surface)
27C soll(klon)--input-R-  Net LW flux above surface (W/m2), positive up
28C                   (i.e. flux emited - flux absorbed by the surface)
29C sens(klon)--input-R-  Sensible Flux at surface  (W/m2), positive down
30C evap(klon)--input-R-  Evaporation + sublimation watter vapour mass flux
31C                   (kg/m2/s), positive up
32C rain_fall(klon)
33C           --input-R- Liquid  watter mass flux (kg/m2/s), positive down
34C snow_fall(klon)
35C           --input-R- Solid  watter mass flux (kg/m2/s), positive down
36C ts(klon)----input-R- Surface temperature (K)
37C d_etp_tot---input-R- Heat flux equivalent to atmospheric enthalpy
38C                    change (W/m2)
39C d_qt_tot----input-R- Mass flux equivalent to atmospheric watter mass
40C                    change (kg/m2/s)
41C d_ec_tot----input-R- Flux equivalent to atmospheric cinetic energy
42C                    change (W/m2)
43C
44C fs_bound---output-R- Thermal flux at the atmosphere boundaries (W/m2)
45C fq_bound---output-R- Watter mass flux at the atmosphere boundaries (kg/m2/s)
46C
47C J.L. Dufresne, July 2002
48C======================================================================
49C
50      use dimphy
51      implicit none
52
53#include "YOMCST.h"
54C
55C     Input variables
56      real airephy(klon)
57      CHARACTER*15 tit
58      INTEGER iprt
59      real tops(klon),topl(klon),sols(klon),soll(klon)
60      real sens(klon),evap(klon),rain_fall(klon),snow_fall(klon)
61      REAL ts(klon)
62      REAL d_etp_tot, d_qt_tot, d_ec_tot
63c     Output variables
64      REAL fs_bound, fq_bound
65C
66C     Local variables
67      real stops,stopl,ssols,ssoll
68      real ssens,sfront,slat
69      real airetot, zcpvap, zcwat, zcice
70      REAL rain_fall_tot, snow_fall_tot, evap_tot
71C
72      integer i
73C
74      integer pas
75      save pas
76      data pas/0/
77C
78      pas=pas+1
79      stops=0.
80      stopl=0.
81      ssols=0.
82      ssoll=0.
83      ssens=0.
84      sfront = 0.
85      evap_tot = 0.
86      rain_fall_tot = 0.
87      snow_fall_tot = 0.
88      airetot=0.
89C
90C     Pour les chaleur specifiques de la vapeur d'eau, de l'eau et de
91C     la glace, on travaille par difference a la chaleur specifique de l'
92c     air sec. En effet, comme on travaille a niveau de pression donne,
93C     toute variation de la masse d'un constituant est totalement
94c     compense par une variation de masse d'air.
95C
96      zcpvap=RCPV-RCPD
97      zcwat=RCW-RCPD
98      zcice=RCS-RCPD
99C
100      do i=1,klon
101           stops=stops+tops(i)*airephy(i)
102           stopl=stopl+topl(i)*airephy(i)
103           ssols=ssols+sols(i)*airephy(i)
104           ssoll=ssoll+soll(i)*airephy(i)
105           ssens=ssens+sens(i)*airephy(i)
106           sfront = sfront
107     $         + ( evap(i)*zcpvap-rain_fall(i)*zcwat-snow_fall(i)*zcice
108     $           ) *ts(i) *airephy(i)
109           evap_tot = evap_tot + evap(i)*airephy(i)
110           rain_fall_tot = rain_fall_tot + rain_fall(i)*airephy(i)
111           snow_fall_tot = snow_fall_tot + snow_fall(i)*airephy(i)
112           airetot=airetot+airephy(i)
113      enddo
114      stops=stops/airetot
115      stopl=stopl/airetot
116      ssols=ssols/airetot
117      ssoll=ssoll/airetot
118      ssens=ssens/airetot
119      sfront = sfront/airetot
120      evap_tot = evap_tot /airetot
121      rain_fall_tot = rain_fall_tot/airetot
122      snow_fall_tot = snow_fall_tot/airetot
123C
124      slat = RLVTT * rain_fall_tot + RLSTT * snow_fall_tot
125C     Heat flux at atm. boundaries
126      fs_bound = stops-stopl - (ssols+ssoll)+ssens+sfront
127     $    + slat
128C     Watter flux at atm. boundaries
129      fq_bound = evap_tot - rain_fall_tot -snow_fall_tot
130C
131      IF (iprt.ge.1) write(6,6666)
132     $    tit, pas, fs_bound, d_etp_tot, fq_bound, d_qt_tot
133C
134      IF (iprt.ge.1) write(6,6668)
135     $    tit, pas, d_etp_tot+d_ec_tot-fs_bound, d_qt_tot-fq_bound
136C
137      IF (iprt.ge.2) write(6,6667)
138     $    tit, pas, stops,stopl,ssols,ssoll,ssens,slat,evap_tot
139     $    ,rain_fall_tot+snow_fall_tot
140
141      return
142
143 6666 format('Phys. Flux Budget ',a15,1i6,x,2(f10.2,x),2(1pE13.5))
144 6667 format('Phys. Boundary Flux ',a15,1i6,x,6(f10.2,x),2(1pE13.5))
145 6668 format('Phys. Total Budget ',a15,1i6,x,f10.2,2(1pE13.5))
146
147      end
148
149C======================================================================
150      SUBROUTINE diagetpq(airephy,tit,iprt,idiag,idiag2,dtime
151     e  ,t,q,ql,qs,u,v,paprs,pplay
152     s  , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
153C======================================================================
154C
155C Purpose:
156C    Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
157C    et calcul le flux de chaleur et le flux d'eau necessaire a ces
158C    changements. Ces valeurs sont moyennees sur la surface de tout
159C    le globe et sont exprime en W/2 et kg/s/m2
160C    Outil pour diagnostiquer la conservation de l'energie
161C    et de la masse dans la physique. Suppose que les niveau de
162c    pression entre couche ne varie pas entre 2 appels.
163C
164C Plusieurs de ces diagnostics peuvent etre fait en parallele: les
165c bilans sont sauvegardes dans des tableaux indices. On parlera
166C "d'indice de diagnostic"
167c
168C
169c======================================================================
170C Arguments:
171C airephy-------input-R-  grid area
172C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
173C iprt----input-I-  PRINT level ( <=1 : no PRINT)
174C idiag---input-I- indice dans lequel sera range les nouveaux
175C                  bilans d' entalpie et de masse
176C idiag2--input-I-les nouveaux bilans d'entalpie et de masse
177C                 sont compare au bilan de d'enthalpie de masse de
178C                 l'indice numero idiag2
179C                 Cas parriculier : si idiag2=0, pas de comparaison, on
180c                 sort directement les bilans d'enthalpie et de masse
181C dtime----input-R- time step (s)
182c t--------input-R- temperature (K)
183c q--------input-R- vapeur d'eau (kg/kg)
184c ql-------input-R- liquid watter (kg/kg)
185c qs-------input-R- solid watter (kg/kg)
186c u--------input-R- vitesse u
187c v--------input-R- vitesse v
188c paprs----input-R- pression a intercouche (Pa)
189c pplay----input-R- pression au milieu de couche (Pa)
190c
191C the following total value are computed by UNIT of earth surface
192C
193C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
194c            change (J/m2) during one time step (dtime) for the whole
195C            atmosphere (air, watter vapour, liquid and solid)
196C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
197C           total watter (kg/m2) change during one time step (dtime),
198C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
199C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
200C d_qs------output-R- same, for the solid watter only (kg/m2/s)
201C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
202C
203C     other (COMMON...)
204C     RCPD, RCPV, ....
205C
206C J.L. Dufresne, July 2002
207c======================================================================
208 
209      use dimphy
210      use cpdet_phy_mod, only: cpdet
211      IMPLICIT NONE
212C
213#include "YOMCST.h"
214C
215c     Input variables
216      real airephy(klon)
217      CHARACTER*15 tit
218      INTEGER iprt,idiag, idiag2
219      REAL dtime
220      REAL t(klon,klev), q(klon,klev), ql(klon,klev), qs(klon,klev)
221      REAL u(klon,klev), v(klon,klev)
222      REAL paprs(klon,klev+1), pplay(klon,klev)
223c     Output variables
224      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
225C
226C     Local variables
227c
228      REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
229     .  , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
230c h_vcol_tot--  total enthalpy of vertical air column
231C            (air with watter vapour, liquid and solid) (J/m2)
232c h_dair_tot-- total enthalpy of dry air (J/m2)
233c h_qw_tot----  total enthalpy of watter vapour (J/m2)
234c h_ql_tot----  total enthalpy of liquid watter (J/m2)
235c h_qs_tot----  total enthalpy of solid watter  (J/m2)
236c qw_tot------  total mass of watter vapour (kg/m2)
237c ql_tot------  total mass of liquid watter (kg/m2)
238c qs_tot------  total mass of solid watter (kg/m2)
239c ec_tot------  total cinetic energy (kg/m2)
240C
241      REAL zairm(klon,klev) ! layer air mass (kg/m2)
242      REAL  zqw_col(klon)
243      REAL  zql_col(klon)
244      REAL  zqs_col(klon)
245      REAL  zec_col(klon)
246      REAL  zh_dair_col(klon)
247      REAL  zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
248C
249      REAL      d_h_dair, d_h_qw, d_h_ql, d_h_qs
250C
251      REAL airetot, zcpvap, zcwat, zcice
252C
253      INTEGER i, k
254C
255      INTEGER ndiag     ! max number of diagnostic in parallel
256      PARAMETER (ndiag=10)
257      integer pas(ndiag)
258      save pas
259      data pas/ndiag*0/
260C     
261      REAL      h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
262     $    , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
263     $    , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
264      SAVE      h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
265     $        , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
266
267c======================================================================
268C
269      DO k = 1, klev
270        DO i = 1, klon
271C         layer air mass
272          zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
273        ENDDO
274      END DO
275C
276C     Reset variables
277      DO i = 1, klon
278        zqw_col(i)=0.
279        zql_col(i)=0.
280        zqs_col(i)=0.
281        zec_col(i) = 0.
282        zh_dair_col(i) = 0.
283        zh_qw_col(i) = 0.
284        zh_ql_col(i) = 0.
285        zh_qs_col(i) = 0.
286      ENDDO
287C
288      zcpvap=RCPV
289      zcwat=RCW
290      zcice=RCS
291C
292C     Compute vertical sum for each atmospheric column
293C     ================================================
294      DO k = 1, klev
295        DO i = 1, klon
296C         Watter mass
297          zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)
298          zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)
299          zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)
300C         Cinetic Energy
301          zec_col(i) =  zec_col(i)
302     $        +0.5*(u(i,k)**2+v(i,k)**2)*zairm(i,k)
303C         Air enthalpy
304! ADAPTATION GCM POUR CP(T)
305          zh_dair_col(i) = zh_dair_col(i)
306     $    + cpdet(t(i,k))*(1.-q(i,k)-ql(i,k)-qs(i,k))*zairm(i,k)*t(i,k)
307          zh_qw_col(i) = zh_qw_col(i)
308     $        + zcpvap*q(i,k)*zairm(i,k)*t(i,k)
309          zh_ql_col(i) = zh_ql_col(i)
310     $        + zcwat*ql(i,k)*zairm(i,k)*t(i,k)
311     $        - RLVTT*ql(i,k)*zairm(i,k)
312          zh_qs_col(i) = zh_qs_col(i)
313     $        + zcice*qs(i,k)*zairm(i,k)*t(i,k)
314     $        - RLSTT*qs(i,k)*zairm(i,k)
315        END DO
316      ENDDO
317C
318C     Mean over the planete surface
319C     =============================
320      qw_tot = 0.
321      ql_tot = 0.
322      qs_tot = 0.
323      ec_tot = 0.
324      h_vcol_tot = 0.
325      h_dair_tot = 0.
326      h_qw_tot = 0.
327      h_ql_tot = 0.
328      h_qs_tot = 0.
329      airetot=0.
330C
331      do i=1,klon
332        qw_tot = qw_tot + zqw_col(i)*airephy(i)
333        ql_tot = ql_tot + zql_col(i)*airephy(i)
334        qs_tot = qs_tot + zqs_col(i)*airephy(i)
335        ec_tot = ec_tot + zec_col(i)*airephy(i)
336        h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
337        h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
338        h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
339        h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
340        airetot=airetot+airephy(i)
341      END DO
342C
343      qw_tot = qw_tot/airetot
344      ql_tot = ql_tot/airetot
345      qs_tot = qs_tot/airetot
346      ec_tot = ec_tot/airetot
347      h_dair_tot = h_dair_tot/airetot
348      h_qw_tot = h_qw_tot/airetot
349      h_ql_tot = h_ql_tot/airetot
350      h_qs_tot = h_qs_tot/airetot
351C
352      h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
353c     print*,'airetot=',airetot,'   h_dair_tot=',h_dair_tot
354C
355C     Compute the change of the atmospheric state compare to the one
356C     stored in "idiag2", and convert it in flux. THis computation
357C     is performed IF idiag2 /= 0 and IF it is not the first CALL
358c     for "idiag"
359C     ===================================
360C
361      IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
362        d_h_vcol  = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
363        d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
364        d_h_qw   = (h_qw_tot  - h_qw_pre(idiag2)  )/dtime
365        d_h_ql   = (h_ql_tot  - h_ql_pre(idiag2)  )/dtime
366        d_h_qs   = (h_qs_tot  - h_qs_pre(idiag2)  )/dtime
367        d_qw     = (qw_tot    - qw_pre(idiag2)    )/dtime
368        d_ql     = (ql_tot    - ql_pre(idiag2)    )/dtime
369        d_qs     = (qs_tot    - qs_pre(idiag2)    )/dtime
370        d_ec     = (ec_tot    - ec_pre(idiag2)    )/dtime
371        d_qt = d_qw + d_ql + d_qs
372      ELSE
373        d_h_vcol = 0.
374        d_h_dair = 0.
375        d_h_qw   = 0.
376        d_h_ql   = 0.
377        d_h_qs   = 0.
378        d_qw     = 0.
379        d_ql     = 0.
380        d_qs     = 0.
381        d_ec     = 0.
382        d_qt     = 0.
383      ENDIF
384C
385      IF (iprt.ge.2) THEN
386        WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
387 9000   format('Phys. Watter Mass Budget (kg/m2/s)',A15
388     $      ,1i6,10(1pE14.6))
389        WRITE(6,9001) tit,pas(idiag), d_h_vcol, h_vcol_tot/dtime
390 9001   format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(E14.6,x))
391        WRITE(6,9002) tit,pas(idiag), d_ec
392 9002   format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F10.2))
393      END IF
394C
395C     Store the new atmospheric state in "idiag"
396C
397      pas(idiag)=pas(idiag)+1
398      h_vcol_pre(idiag)  = h_vcol_tot
399      h_dair_pre(idiag) = h_dair_tot
400      h_qw_pre(idiag)   = h_qw_tot
401      h_ql_pre(idiag)   = h_ql_tot
402      h_qs_pre(idiag)   = h_qs_tot
403      qw_pre(idiag)     = qw_tot
404      ql_pre(idiag)     = ql_tot
405      qs_pre(idiag)     = qs_tot
406      ec_pre (idiag)    = ec_tot
407C
408      RETURN
409      END
Note: See TracBrowser for help on using the repository browser.