trunk/phylmd/diagetpq.f90

module diagetpq_m

  IMPLICIT NONE

contains

  SUBROUTINE diagetpq(airephy, tit, iprt, idiag, idiag2, dtime, t, q, ql, qs, &
       u, v, paprs, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)

    ! From LMDZ4/libf/phylmd/diagphy.F, version 1.1.1.1 2004/05/19 12:53:08

    ! Purpose:

    ! Calcule la différence d'enthalpie et de masse d'eau entre deux
    ! appels et calcule le flux de chaleur et le flux d'eau
    ! nécessaires à ces changements. Ces valeurs sont moyennées sur la
    ! surface de tout le globe et sont exprimées en W/m2 et
    ! kg/s/m2. Outil pour diagnostiquer la conservation de l'énergie
    ! et de la masse dans la physique. Suppose que les niveaux de
    ! pression entre les couches ne varient pas entre deux appels.

    ! Plusieurs de ces diagnostics peuvent être faits en parallèle :
    ! les bilans sont sauvegardés dans des tableaux indices. On
    ! parlera "d'indice de diagnostic".

    ! Jean-Louis Dufresne, July 2002

    USE dimphy, ONLY: klev, klon
    USE suphec_m, ONLY: rcpd, rcpv, rcs, rcw, rg, rlstt, rlvtt

    ! Arguments: 
    ! airephy-------input-R- grid area
    ! iprt----input-I- PRINT level ( <=1 : no PRINT)
    ! idiag---input-I- indice dans lequel sera range les nouveaux
    ! bilans d' entalpie et de masse
    ! idiag2--input-I-les nouveaux bilans d'entalpie et de masse 
    ! sont compare au bilan de d'enthalpie de masse de
    ! l'indice numero idiag2 
    ! Cas particulier : si idiag2=0, pas de comparaison, on
    ! sort directement les bilans d'enthalpie et de masse 
    ! dtime----input-R- time step (s)
    ! t--------input-R- temperature (K)
    ! q--------input-R- vapeur d'eau (kg/kg)
    ! ql-------input-R- liquid water (kg/kg)
    ! qs-------input-R- solid water (kg/kg)
    ! u--------input-R- vitesse u
    ! v--------input-R- vitesse v
    ! paprs----input-R- pression a intercouche (Pa)

    ! the following total value are computed by UNIT of earth surface

    ! d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy 
    ! change (J/m2) during one time step (dtime) for the whole 
    ! atmosphere (air, water vapour, liquid and solid)
    ! d_qt------output-R- total water mass flux (kg/m2/s) defined as the 
    ! total water (kg/m2) change during one time step (dtime),
    ! d_qw------output-R- same, for the water vapour only (kg/m2/s)
    ! d_ql------output-R- same, for the liquid water only (kg/m2/s)
    ! d_qs------output-R- same, for the solid water only (kg/m2/s)
    ! d_ec------output-R- Kinetic Energy Budget (W/m2) for vertical air column

    ! other (COMMON...)
    ! RCPD, RCPV, ....

    ! Input variables
    real airephy(klon)
    CHARACTER(len=*), intent(in):: tit ! comment added in PRINT
    INTEGER iprt, idiag, idiag2
    REAL, intent(in):: dtime
    REAL, intent(in):: t(klon, klev)
    REAL, intent(in):: q(klon, klev), ql(klon, klev), qs(klon, klev)
    REAL u(klon, klev), v(klon, klev)
    REAL, intent(in):: paprs(klon, klev+1)
    ! Output variables
    REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec

    ! Local variables

    REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot, h_qs_tot, qw_tot, ql_tot
    real qs_tot , ec_tot
    ! h_vcol_tot-- total enthalpy of vertical air column 
    ! (air with water vapour, liquid and solid) (J/m2)
    ! h_dair_tot-- total enthalpy of dry air (J/m2)
    ! h_qw_tot---- total enthalpy of water vapour (J/m2)
    ! h_ql_tot---- total enthalpy of liquid water (J/m2)
    ! h_qs_tot---- total enthalpy of solid water (J/m2)
    ! qw_tot------ total mass of water vapour (kg/m2)
    ! ql_tot------ total mass of liquid water (kg/m2)
    ! qs_tot------ total mass of solid water (kg/m2)
    ! ec_tot------ total kinetic energy (kg/m2)

    REAL zairm(klon, klev) ! layer air mass (kg/m2)
    REAL zqw_col(klon)
    REAL zql_col(klon)
    REAL zqs_col(klon)
    REAL zec_col(klon)
    REAL zh_dair_col(klon)
    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

    REAL airetot, zcpvap, zcwat, zcice

    INTEGER i, k

    INTEGER, PARAMETER:: ndiag = 10 ! max number of diagnostic in parallel
    integer:: pas(ndiag) = 0

    REAL, save:: h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
    REAL, save:: h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag), ql_pre(ndiag)
    REAL, save:: qs_pre(ndiag), ec_pre(ndiag)

    !-------------------------------------------------------------

    DO k = 1, klev
       DO i = 1, klon
          ! layer air mass
          zairm(i, k) = (paprs(i, k)-paprs(i, k+1))/RG
       ENDDO
    END DO

    ! Reset variables
    DO i = 1, klon
       zqw_col(i)=0.
       zql_col(i)=0.
       zqs_col(i)=0.
       zec_col(i) = 0.
       zh_dair_col(i) = 0.
       zh_qw_col(i) = 0.
       zh_ql_col(i) = 0.
       zh_qs_col(i) = 0.
    ENDDO

    zcpvap=RCPV
    zcwat=RCW
    zcice=RCS

    ! Compute vertical sum for each atmospheric column
    DO k = 1, klev
       DO i = 1, klon
          ! Water mass
          zqw_col(i) = zqw_col(i) + q(i, k)*zairm(i, k)
          zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
          zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
          ! Kinetic Energy
          zec_col(i) = zec_col(i) +0.5*(u(i, k)**2+v(i, k)**2)*zairm(i, k)
          ! Air enthalpy
          zh_dair_col(i) = zh_dair_col(i) &
               + RCPD*(1.-q(i, k)-ql(i, k)-qs(i, k))*zairm(i, k)*t(i, k)
          zh_qw_col(i) = zh_qw_col(i) + zcpvap*q(i, k)*zairm(i, k)*t(i, k) 
          zh_ql_col(i) = zh_ql_col(i) &
               + zcwat*ql(i, k)*zairm(i, k)*t(i, k) &
               - RLVTT*ql(i, k)*zairm(i, k)
          zh_qs_col(i) = zh_qs_col(i) &
               + zcice*qs(i, k)*zairm(i, k)*t(i, k) &
               - RLSTT*qs(i, k)*zairm(i, k)
       END DO
    ENDDO

    ! Mean over the planet surface
    qw_tot = 0.
    ql_tot = 0.
    qs_tot = 0.
    ec_tot = 0.
    h_vcol_tot = 0.
    h_dair_tot = 0.
    h_qw_tot = 0.
    h_ql_tot = 0.
    h_qs_tot = 0.
    airetot=0.

    do i=1, klon
       qw_tot = qw_tot + zqw_col(i)*airephy(i)
       ql_tot = ql_tot + zql_col(i)*airephy(i)
       qs_tot = qs_tot + zqs_col(i)*airephy(i)
       ec_tot = ec_tot + zec_col(i)*airephy(i)
       h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
       h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
       h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
       h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
       airetot=airetot+airephy(i)
    END DO

    qw_tot = qw_tot/airetot
    ql_tot = ql_tot/airetot
    qs_tot = qs_tot/airetot
    ec_tot = ec_tot/airetot
    h_dair_tot = h_dair_tot/airetot
    h_qw_tot = h_qw_tot/airetot
    h_ql_tot = h_ql_tot/airetot
    h_qs_tot = h_qs_tot/airetot

    h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot

    ! Compute the change of the atmospheric state compared to the one
    ! stored in "idiag2", and convert it in flux. This computation is
    ! performed if idiag2 /= 0 and if it is not the first call for
    ! "idiag".

    IF ((idiag2 > 0) .and. (pas(idiag2) /= 0)) THEN
       d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
       d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
       d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
       d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime 
       d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime 
       d_qw = (qw_tot - qw_pre(idiag2) )/dtime
       d_ql = (ql_tot - ql_pre(idiag2) )/dtime
       d_qs = (qs_tot - qs_pre(idiag2) )/dtime
       d_ec = (ec_tot - ec_pre(idiag2) )/dtime
       d_qt = d_qw + d_ql + d_qs
    ELSE 
       d_h_vcol = 0.
       d_h_dair = 0.
       d_h_qw = 0.
       d_h_ql = 0.
       d_h_qs = 0. 
       d_qw = 0.
       d_ql = 0.
       d_qs = 0.
       d_ec = 0.
       d_qt = 0.
    ENDIF

    IF (iprt >= 2) THEN
       WRITE(6, 9000) tit, pas(idiag), d_qt, d_qw, d_ql, d_qs
9000   format('Phys. Water Mass Budget (kg/m2/s)', A15, 1i6, 10(1pE14.6))
       WRITE(6, 9001) tit, pas(idiag), d_h_vcol
9001   format('Phys. Enthalpy Budget (W/m2) ', A15, 1i6, 10(F8.2))
       WRITE(6, 9002) tit, pas(idiag), d_ec
9002   format('Phys. Kinetic Energy Budget (W/m2) ', A15, 1i6, 10(F8.2))
    END IF

    ! Store the new atmospheric state in "idiag"
    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

  END SUBROUTINE diagetpq

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