/[lmdze]/trunk/libf/phylmd/diagetpq.f90
ViewVC logotype

Diff of /trunk/libf/phylmd/diagetpq.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/diagetpq.f revision 47 by guez, Fri Jul 1 15:00:48 2011 UTC trunk/libf/phylmd/diagetpq.f90 revision 52 by guez, Fri Sep 23 12:28:01 2011 UTC
# Line 1  Line 1 
1        SUBROUTINE diagetpq(airephy,tit,iprt,idiag,idiag2,dtime  module diagetpq_m
2       e  ,t,q,ql,qs,u,v,paprs  
3       s  , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)    IMPLICIT NONE
4  !  
5  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/diagphy.F,v 1.1.1.1 2004/05/19 12:53:08 lmdzadmin Exp $  contains
6  !  
7  C======================================================================    SUBROUTINE diagetpq(airephy, tit, iprt, idiag, idiag2, dtime, t, q, ql, qs, &
8  C         u, v, paprs, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
9  C Purpose:  
10  C    Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,      ! From LMDZ4/libf/phylmd/diagphy.F, version 1.1.1.1 2004/05/19 12:53:08
11  C    et calcul le flux de chaleur et le flux d'eau necessaire a ces  
12  C    changements. Ces valeurs sont moyennees sur la surface de tout      ! Purpose:
13  C    le globe et sont exprime en W/2 et kg/s/m2  
14  C    Outil pour diagnostiquer la conservation de l'energie      ! Calcule la différence d'enthalpie et de masse d'eau entre deux
15  C    et de la masse dans la physique. Suppose que les niveau de      ! appels et calcule le flux de chaleur et le flux d'eau
16  c    pression entre couche ne varie pas entre 2 appels.      ! nécessaires à ces changements. Ces valeurs sont moyennées sur la
17  C      ! surface de tout le globe et sont exprimées en W/m2 et
18  C Plusieurs de ces diagnostics peuvent etre fait en parallele: les      ! kg/s/m2. Outil pour diagnostiquer la conservation de l'énergie
19  c bilans sont sauvegardes dans des tableaux indices. On parlera      ! et de la masse dans la physique. Suppose que les niveaux de
20  C "d'indice de diagnostic"      ! pression entre les couches ne varient pas entre deux appels.
21  c  
22  C      ! Plusieurs de ces diagnostics peuvent être faits en parallèle :
23  c======================================================================      ! les bilans sont sauvegardés dans des tableaux indices. On
24  C Arguments:      ! parlera "d'indice de diagnostic".
25  C airephy-------input-R-  grid area  
26  C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)      ! Jean-Louis Dufresne, July 2002
27  C iprt----input-I-  PRINT level ( <=1 : no PRINT)  
28  C idiag---input-I- indice dans lequel sera range les nouveaux      USE dimphy, ONLY: klev, klon
29  C                  bilans d' entalpie et de masse      USE suphec_m, ONLY: rcpd, rcpv, rcs, rcw, rg, rlstt, rlvtt
30  C idiag2--input-I-les nouveaux bilans d'entalpie et de masse  
31  C                 sont compare au bilan de d'enthalpie de masse de      ! Arguments:
32  C                 l'indice numero idiag2      ! airephy-------input-R- grid area
33  C                 Cas parriculier : si idiag2=0, pas de comparaison, on      ! tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
34  c                 sort directement les bilans d'enthalpie et de masse      ! iprt----input-I- PRINT level ( <=1 : no PRINT)
35  C dtime----input-R- time step (s)      ! idiag---input-I- indice dans lequel sera range les nouveaux
36  c t--------input-R- temperature (K)      ! bilans d' entalpie et de masse
37  c q--------input-R- vapeur d'eau (kg/kg)      ! idiag2--input-I-les nouveaux bilans d'entalpie et de masse
38  c ql-------input-R- liquid watter (kg/kg)      ! sont compare au bilan de d'enthalpie de masse de
39  c qs-------input-R- solid watter (kg/kg)      ! l'indice numero idiag2
40  c u--------input-R- vitesse u      ! Cas particulier : si idiag2=0, pas de comparaison, on
41  c v--------input-R- vitesse v      ! sort directement les bilans d'enthalpie et de masse
42  c paprs----input-R- pression a intercouche (Pa)      ! dtime----input-R- time step (s)
43  c      ! t--------input-R- temperature (K)
44  C the following total value are computed by UNIT of earth surface      ! q--------input-R- vapeur d'eau (kg/kg)
45  C      ! ql-------input-R- liquid water (kg/kg)
46  C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy      ! qs-------input-R- solid water (kg/kg)
47  c            change (J/m2) during one time step (dtime) for the whole      ! u--------input-R- vitesse u
48  C            atmosphere (air, watter vapour, liquid and solid)      ! v--------input-R- vitesse v
49  C d_qt------output-R- total water mass flux (kg/m2/s) defined as the      ! paprs----input-R- pression a intercouche (Pa)
50  C           total watter (kg/m2) change during one time step (dtime),  
51  C d_qw------output-R- same, for the watter vapour only (kg/m2/s)      ! the following total value are computed by UNIT of earth surface
52  C d_ql------output-R- same, for the liquid watter only (kg/m2/s)  
53  C d_qs------output-R- same, for the solid watter only (kg/m2/s)      ! d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
54  C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column      ! change (J/m2) during one time step (dtime) for the whole
55  C      ! atmosphere (air, water vapour, liquid and solid)
56  C     other (COMMON...)      ! d_qt------output-R- total water mass flux (kg/m2/s) defined as the
57  C     RCPD, RCPV, ....      ! total water (kg/m2) change during one time step (dtime),
58  C      ! d_qw------output-R- same, for the water vapour only (kg/m2/s)
59  C J.L. Dufresne, July 2002      ! d_ql------output-R- same, for the liquid water only (kg/m2/s)
60  c======================================================================      ! d_qs------output-R- same, for the solid water only (kg/m2/s)
61        ! d_ec------output-R- Kinetic Energy Budget (W/m2) for vertical air column
62        use dimens_m  
63        use dimphy      ! other (COMMON...)
64        use SUPHEC_M      ! RCPD, RCPV, ....
65        use yoethf_m  
66        IMPLICIT NONE      ! Input variables
67  C      real airephy(klon)
68  C      CHARACTER(len=15) tit
69  c     Input variables      INTEGER iprt, idiag, idiag2
70        real airephy(klon)      REAL, intent(in):: dtime
71        CHARACTER*15 tit      REAL, intent(in):: t(klon, klev)
72        INTEGER iprt,idiag, idiag2      REAL, intent(in):: q(klon, klev), ql(klon, klev), qs(klon, klev)
73        REAL, intent(in):: dtime      REAL u(klon, klev), v(klon, klev)
74        REAL t(klon,klev), q(klon,klev), ql(klon,klev), qs(klon,klev)      REAL, intent(in):: paprs(klon, klev+1)
75        REAL u(klon,klev), v(klon,klev)      ! Output variables
76        REAL, intent(in):: paprs(klon,klev+1)      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
77  c     Output variables  
78        REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      ! Local variables
79  C  
80  C     Local variables      REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot, h_qs_tot, qw_tot, ql_tot
81  c      real qs_tot , ec_tot
82        REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot      ! h_vcol_tot-- total enthalpy of vertical air column
83       .  , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot      ! (air with water vapour, liquid and solid) (J/m2)
84  c h_vcol_tot--  total enthalpy of vertical air column      ! h_dair_tot-- total enthalpy of dry air (J/m2)
85  C            (air with watter vapour, liquid and solid) (J/m2)      ! h_qw_tot---- total enthalpy of water vapour (J/m2)
86  c h_dair_tot-- total enthalpy of dry air (J/m2)      ! h_ql_tot---- total enthalpy of liquid water (J/m2)
87  c h_qw_tot----  total enthalpy of watter vapour (J/m2)      ! h_qs_tot---- total enthalpy of solid water (J/m2)
88  c h_ql_tot----  total enthalpy of liquid watter (J/m2)      ! qw_tot------ total mass of water vapour (kg/m2)
89  c h_qs_tot----  total enthalpy of solid watter  (J/m2)      ! ql_tot------ total mass of liquid water (kg/m2)
90  c qw_tot------  total mass of watter vapour (kg/m2)      ! qs_tot------ total mass of solid water (kg/m2)
91  c ql_tot------  total mass of liquid watter (kg/m2)      ! ec_tot------ total kinetic energy (kg/m2)
92  c qs_tot------  total mass of solid watter (kg/m2)  
93  c ec_tot------  total cinetic energy (kg/m2)      REAL zairm(klon, klev) ! layer air mass (kg/m2)
94  C      REAL zqw_col(klon)
95        REAL zairm(klon,klev) ! layer air mass (kg/m2)      REAL zql_col(klon)
96        REAL  zqw_col(klon)      REAL zqs_col(klon)
97        REAL  zql_col(klon)      REAL zec_col(klon)
98        REAL  zqs_col(klon)      REAL zh_dair_col(klon)
99        REAL  zec_col(klon)      REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
100        REAL  zh_dair_col(klon)  
101        REAL  zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)      REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
102  C  
103        REAL      d_h_dair, d_h_qw, d_h_ql, d_h_qs      REAL airetot, zcpvap, zcwat, zcice
104  C  
105        REAL airetot, zcpvap, zcwat, zcice      INTEGER i, k
106  C  
107        INTEGER i, k      INTEGER, PARAMETER:: ndiag = 10 ! max number of diagnostic in parallel
108  C      integer:: pas(ndiag) = 0
109        INTEGER ndiag     ! max number of diagnostic in parallel  
110        PARAMETER (ndiag=10)      REAL, save:: h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
111        integer pas(ndiag)      REAL, save:: h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag), ql_pre(ndiag)
112        save pas      REAL, save:: qs_pre(ndiag), ec_pre(ndiag)
113        data pas/ndiag*0/  
114  C          !-------------------------------------------------------------
115        REAL      h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)  
116       $    , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)      DO k = 1, klev
117       $    , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)         DO i = 1, klon
118        SAVE      h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre            ! layer air mass
119       $        , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre            zairm(i, k) = (paprs(i, k)-paprs(i, k+1))/RG
120           ENDDO
121  c======================================================================      END DO
122  C  
123        DO k = 1, klev      ! Reset variables
124          DO i = 1, klon      DO i = 1, klon
125  C         layer air mass         zqw_col(i)=0.
126            zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG         zql_col(i)=0.
127          ENDDO         zqs_col(i)=0.
128        END DO         zec_col(i) = 0.
129  C         zh_dair_col(i) = 0.
130  C     Reset variables         zh_qw_col(i) = 0.
131        DO i = 1, klon         zh_ql_col(i) = 0.
132          zqw_col(i)=0.         zh_qs_col(i) = 0.
133          zql_col(i)=0.      ENDDO
134          zqs_col(i)=0.  
135          zec_col(i) = 0.      zcpvap=RCPV
136          zh_dair_col(i) = 0.      zcwat=RCW
137          zh_qw_col(i) = 0.      zcice=RCS
138          zh_ql_col(i) = 0.  
139          zh_qs_col(i) = 0.      ! Compute vertical sum for each atmospheric column
140        ENDDO      DO k = 1, klev
141  C         DO i = 1, klon
142        zcpvap=RCPV            ! Water mass
143        zcwat=RCW            zqw_col(i) = zqw_col(i) + q(i, k)*zairm(i, k)
144        zcice=RCS            zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
145  C            zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
146  C     Compute vertical sum for each atmospheric column            ! Kinetic Energy
147  C     ================================================            zec_col(i) = zec_col(i) +0.5*(u(i, k)**2+v(i, k)**2)*zairm(i, k)
148        DO k = 1, klev            ! Air enthalpy
149          DO i = 1, klon            zh_dair_col(i) = zh_dair_col(i) &
150  C         Watter mass                 + RCPD*(1.-q(i, k)-ql(i, k)-qs(i, k))*zairm(i, k)*t(i, k)
151            zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)            zh_qw_col(i) = zh_qw_col(i) + zcpvap*q(i, k)*zairm(i, k)*t(i, k)
152            zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)            zh_ql_col(i) = zh_ql_col(i) &
153            zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)                 + zcwat*ql(i, k)*zairm(i, k)*t(i, k) &
154  C         Cinetic Energy                 - RLVTT*ql(i, k)*zairm(i, k)
155            zec_col(i) =  zec_col(i)            zh_qs_col(i) = zh_qs_col(i) &
156       $        +0.5*(u(i,k)**2+v(i,k)**2)*zairm(i,k)                 + zcice*qs(i, k)*zairm(i, k)*t(i, k) &
157  C         Air enthalpy                 - RLSTT*qs(i, k)*zairm(i, k)
158            zh_dair_col(i) = zh_dair_col(i)         END DO
159       $        + RCPD*(1.-q(i,k)-ql(i,k)-qs(i,k))*zairm(i,k)*t(i,k)      ENDDO
160            zh_qw_col(i) = zh_qw_col(i)  
161       $        + zcpvap*q(i,k)*zairm(i,k)*t(i,k)      ! Mean over the planet surface
162            zh_ql_col(i) = zh_ql_col(i)      qw_tot = 0.
163       $        + zcwat*ql(i,k)*zairm(i,k)*t(i,k)      ql_tot = 0.
164       $        - RLVTT*ql(i,k)*zairm(i,k)      qs_tot = 0.
165            zh_qs_col(i) = zh_qs_col(i)      ec_tot = 0.
166       $        + zcice*qs(i,k)*zairm(i,k)*t(i,k)      h_vcol_tot = 0.
167       $        - RLSTT*qs(i,k)*zairm(i,k)      h_dair_tot = 0.
168        h_qw_tot = 0.
169          END DO      h_ql_tot = 0.
170        ENDDO      h_qs_tot = 0.
171  C      airetot=0.
172  C     Mean over the planete surface  
173  C     =============================      do i=1, klon
174        qw_tot = 0.         qw_tot = qw_tot + zqw_col(i)*airephy(i)
175        ql_tot = 0.         ql_tot = ql_tot + zql_col(i)*airephy(i)
176        qs_tot = 0.         qs_tot = qs_tot + zqs_col(i)*airephy(i)
177        ec_tot = 0.         ec_tot = ec_tot + zec_col(i)*airephy(i)
178        h_vcol_tot = 0.         h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
179        h_dair_tot = 0.         h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
180        h_qw_tot = 0.         h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
181        h_ql_tot = 0.         h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
182        h_qs_tot = 0.         airetot=airetot+airephy(i)
183        airetot=0.      END DO
184  C  
185        do i=1,klon      qw_tot = qw_tot/airetot
186          qw_tot = qw_tot + zqw_col(i)*airephy(i)      ql_tot = ql_tot/airetot
187          ql_tot = ql_tot + zql_col(i)*airephy(i)      qs_tot = qs_tot/airetot
188          qs_tot = qs_tot + zqs_col(i)*airephy(i)      ec_tot = ec_tot/airetot
189          ec_tot = ec_tot + zec_col(i)*airephy(i)      h_dair_tot = h_dair_tot/airetot
190          h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)      h_qw_tot = h_qw_tot/airetot
191          h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)      h_ql_tot = h_ql_tot/airetot
192          h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)      h_qs_tot = h_qs_tot/airetot
193          h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)  
194          airetot=airetot+airephy(i)      h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
195        END DO  
196  C      ! Compute the change of the atmospheric state compared to the one
197        qw_tot = qw_tot/airetot      ! stored in "idiag2", and convert it in flux. This computation is
198        ql_tot = ql_tot/airetot      ! performed if idiag2 /= 0 and if it is not the first call for
199        qs_tot = qs_tot/airetot      ! "idiag".
200        ec_tot = ec_tot/airetot  
201        h_dair_tot = h_dair_tot/airetot      IF ((idiag2 > 0) .and. (pas(idiag2) /= 0)) THEN
202        h_qw_tot = h_qw_tot/airetot         d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
203        h_ql_tot = h_ql_tot/airetot         d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
204        h_qs_tot = h_qs_tot/airetot         d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
205  C         d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
206        h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot         d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
207  C         d_qw = (qw_tot - qw_pre(idiag2) )/dtime
208  C     Compute the change of the atmospheric state compare to the one         d_ql = (ql_tot - ql_pre(idiag2) )/dtime
209  C     stored in "idiag2", and convert it in flux. THis computation         d_qs = (qs_tot - qs_pre(idiag2) )/dtime
210  C     is performed IF idiag2 /= 0 and IF it is not the first CALL         d_ec = (ec_tot - ec_pre(idiag2) )/dtime
211  c     for "idiag"         d_qt = d_qw + d_ql + d_qs
212  C     ===================================      ELSE
213  C         d_h_vcol = 0.
214        IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN         d_h_dair = 0.
215          d_h_vcol  = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime         d_h_qw = 0.
216          d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime         d_h_ql = 0.
217          d_h_qw   = (h_qw_tot  - h_qw_pre(idiag2)  )/dtime         d_h_qs = 0.
218          d_h_ql   = (h_ql_tot  - h_ql_pre(idiag2)  )/dtime         d_qw = 0.
219          d_h_qs   = (h_qs_tot  - h_qs_pre(idiag2)  )/dtime         d_ql = 0.
220          d_qw     = (qw_tot    - qw_pre(idiag2)    )/dtime         d_qs = 0.
221          d_ql     = (ql_tot    - ql_pre(idiag2)    )/dtime         d_ec = 0.
222          d_qs     = (qs_tot    - qs_pre(idiag2)    )/dtime         d_qt = 0.
223          d_ec     = (ec_tot    - ec_pre(idiag2)    )/dtime      ENDIF
224          d_qt = d_qw + d_ql + d_qs  
225        ELSE      IF (iprt >= 2) THEN
226          d_h_vcol = 0.         WRITE(6, 9000) tit, pas(idiag), d_qt, d_qw, d_ql, d_qs
227          d_h_dair = 0.  9000   format('Phys. Water Mass Budget (kg/m2/s)', A15, 1i6, 10(1pE14.6))
228          d_h_qw   = 0.         WRITE(6, 9001) tit, pas(idiag), d_h_vcol
229          d_h_ql   = 0.  9001   format('Phys. Enthalpy Budget (W/m2) ', A15, 1i6, 10(F8.2))
230          d_h_qs   = 0.         WRITE(6, 9002) tit, pas(idiag), d_ec
231          d_qw     = 0.  9002   format('Phys. Kinetic Energy Budget (W/m2) ', A15, 1i6, 10(F8.2))
232          d_ql     = 0.      END IF
233          d_qs     = 0.  
234          d_ec     = 0.      ! Store the new atmospheric state in "idiag"
235          d_qt     = 0.      pas(idiag)=pas(idiag)+1
236        ENDIF      h_vcol_pre(idiag) = h_vcol_tot
237  C      h_dair_pre(idiag) = h_dair_tot
238        IF (iprt.ge.2) THEN      h_qw_pre(idiag) = h_qw_tot
239          WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs      h_ql_pre(idiag) = h_ql_tot
240   9000   format('Phys. Watter Mass Budget (kg/m2/s)',A15      h_qs_pre(idiag) = h_qs_tot
241       $      ,1i6,10(1pE14.6))      qw_pre(idiag) = qw_tot
242          WRITE(6,9001) tit,pas(idiag), d_h_vcol      ql_pre(idiag) = ql_tot
243   9001   format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))      qs_pre(idiag) = qs_tot
244          WRITE(6,9002) tit,pas(idiag), d_ec      ec_pre (idiag) = ec_tot
245   9002   format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))  
246        END IF    END SUBROUTINE diagetpq
247  C  
248  C     Store the new atmospheric state in "idiag"  end module diagetpq_m
 C  
       pas(idiag)=pas(idiag)+1  
       h_vcol_pre(idiag)  = h_vcol_tot  
       h_dair_pre(idiag) = h_dair_tot  
       h_qw_pre(idiag)   = h_qw_tot  
       h_ql_pre(idiag)   = h_ql_tot  
       h_qs_pre(idiag)   = h_qs_tot  
       qw_pre(idiag)     = qw_tot  
       ql_pre(idiag)     = ql_tot  
       qs_pre(idiag)     = qs_tot  
       ec_pre (idiag)    = ec_tot  
 C  
       RETURN  
       END  

Legend:
Removed from v.47  
changed lines
  Added in v.52

  ViewVC Help
Powered by ViewVC 1.1.21