trunk/phylmd/diagetpq.f

module diagetpq_m

  IMPLICIT NONE

contains

  SUBROUTINE diagetpq(airephy, tit, iprt, idiag, idiag2, dtime, t, q, ql, &
       u, v, paprs, d_h_vcol, d_qt, 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: 

    ! Input variables
    real, intent(in):: airephy(klon) ! grid area
    CHARACTER(len=*), intent(in):: tit ! comment added in PRINT
    INTEGER, intent(in):: iprt ! PRINT level ( <=1 : no PRINT)

    INTEGER, intent(in):: idiag
    ! indice dans lequel seront rangés les nouveaux bilans d'enthalpie et
    ! de masse

    INTEGER, intent(in):: idiag2
    ! Les nouveaux bilans d'enthalpie et de masse sont comparés au
    ! bilan de d'enthalpie de masse de l'indice numéro idiag2. Cas
    ! particulier : si idiag2=0, pas de comparaison, on sort
    ! directement les bilans d'enthalpie et de masse.

    REAL, intent(in):: dtime ! time step (s)
    REAL, intent(in):: t(klon, klev) ! temperature (K)
    REAL, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)
    REAL, intent(in):: ql(klon, klev) ! liquid water (kg/kg)
    REAL, intent(in):: u(klon, klev), v(klon, klev) ! vitesse
    REAL, intent(in):: paprs(klon, klev+1) ! pression a intercouche (Pa)

    ! The following total values are computed per UNIT of Earth surface

    REAL, intent(out):: d_h_vcol
    ! heat flux (W/m2) defined as the enthalpy change (J/m2) during
    ! one time step (dtime) for the whole atmosphere (air, water
    ! vapour, liquid and solid)

    REAL, intent(out):: d_qt
    ! total water mass flux (kg/m2/s) defined as the 
    ! total water (kg/m2) change during one time step (dtime)

    REAL, intent(out):: d_ec
    ! kinetic energy budget (W/m2) for vertical air column

    ! Local:

    REAL d_qw
    ! water vapour mass flux (kg/m2/s) defined as the water vapour
    ! (kg/m2) change during one time step (dtime)

    REAL d_ql ! same, for the liquid water only (kg/m2/s)

    REAL h_vcol_tot
    ! total enthalpy of vertical air column (air with water vapour,
    ! liquid and solid) (J/m2)

    REAL h_dair_tot ! total enthalpy of dry air (J/m2)
    REAL h_qw_tot ! total enthalpy of water vapour (J/m2)
    REAL h_ql_tot ! total enthalpy of liquid water (J/m2)
    REAL qw_tot ! total mass of water vapour (kg/m2)
    REAL ql_tot ! total mass of liquid water (kg/m2)
    real 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 zec_col(klon)
    REAL zh_dair_col(klon)
    REAL zh_qw_col(klon), zh_ql_col(klon)
    REAL d_h_dair, d_h_qw, d_h_ql
    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), qw_pre(ndiag), ql_pre(ndiag)
    REAL, save:: ec_pre(ndiag)

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

    DO k = 1, klev
       DO i = 1, klon
          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.
       zec_col(i) = 0.
       zh_dair_col(i) = 0.
       zh_qw_col(i) = 0.
       zh_ql_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)
          ! 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))*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)
       END DO
    ENDDO

    ! Mean over the planet surface
    qw_tot = 0.
    ql_tot = 0.
    ec_tot = 0.
    h_vcol_tot = 0.
    h_dair_tot = 0.
    h_qw_tot = 0.
    h_ql_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)
       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)
       airetot=airetot+airephy(i)
    END DO

    qw_tot = qw_tot/airetot
    ql_tot = ql_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_vcol_tot = h_dair_tot+h_qw_tot+h_ql_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_qw = (qw_tot - qw_pre(idiag2) )/dtime
       d_ql = (ql_tot - ql_pre(idiag2) )/dtime
       d_ec = (ec_tot - ec_pre(idiag2) )/dtime
       d_qt = d_qw + d_ql
    ELSE 
       d_h_vcol = 0.
       d_h_dair = 0.
       d_h_qw = 0.
       d_h_ql = 0.
       d_qw = 0.
       d_ql = 0.
       d_ec = 0.
       d_qt = 0.
    ENDIF

    IF (iprt >= 2) THEN
       print 9000, tit, pas(idiag), d_qt, d_qw, d_ql
       print 9001, tit, pas(idiag), d_h_vcol
       print 9002, tit, pas(idiag), d_ec
    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
    qw_pre(idiag) = qw_tot
    ql_pre(idiag) = ql_tot
    ec_pre (idiag) = ec_tot

9000 format('Physics water mass budget (kg/m2/s)', A15, 1i6, 10(1pE14.6))
9001 format('Physics enthalpy budget (W/m2) ', A15, 1i6, 10(F8.2))
9002 format('Physics kinetic energy budget (W/m2) ', A15, 1i6, 10(F8.2))

  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	guez	98	SUBROUTINE diagetpq(airephy, tit, iprt, idiag, idiag2, dtime, t, q, ql, &
8			u, v, paprs, d_h_vcol, d_qt, d_ec)
9	guez	51
10	guez	91	! From LMDZ4/libf/phylmd/diagphy.F, version 1.1.1.1, 2004/05/19 12:53:08
11	guez	51
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	guez	91
33			! Input variables
34	guez	98	real, intent(in):: airephy(klon) ! grid area
35	guez	91	CHARACTER(len=*), intent(in):: tit ! comment added in PRINT
36	guez	98	INTEGER, intent(in):: iprt ! PRINT level ( <=1 : no PRINT)
37
38			INTEGER, intent(in):: idiag
39			! indice dans lequel seront rangés les nouveaux bilans d'enthalpie et
40			! de masse
41
42			INTEGER, intent(in):: idiag2
43			! Les nouveaux bilans d'enthalpie et de masse sont comparés au
44			! bilan de d'enthalpie de masse de l'indice numéro idiag2. Cas
45			! particulier : si idiag2=0, pas de comparaison, on sort
46			! directement les bilans d'enthalpie et de masse.
47
48			REAL, intent(in):: dtime ! time step (s)
49			REAL, intent(in):: t(klon, klev) ! temperature (K)
50			REAL, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)
51			REAL, intent(in):: ql(klon, klev) ! liquid water (kg/kg)
52	guez	91	REAL, intent(in):: u(klon, klev), v(klon, klev) ! vitesse
53			REAL, intent(in):: paprs(klon, klev+1) ! pression a intercouche (Pa)
54	guez	51
55	guez	98	! The following total values are computed per UNIT of Earth surface
56	guez	51
57	guez	98	REAL, intent(out):: d_h_vcol
58			! heat flux (W/m2) defined as the enthalpy change (J/m2) during
59			! one time step (dtime) for the whole atmosphere (air, water
60			! vapour, liquid and solid)
61	guez	51
62	guez	98	REAL, intent(out):: d_qt
63			! total water mass flux (kg/m2/s) defined as the
64			! total water (kg/m2) change during one time step (dtime)
65	guez	51
66	guez	98	REAL, intent(out):: d_ec
67			! kinetic energy budget (W/m2) for vertical air column
68	guez	51
69	guez	98	! Local:
70
71			REAL d_qw
72			! water vapour mass flux (kg/m2/s) defined as the water vapour
73			! (kg/m2) change during one time step (dtime)
74
75			REAL d_ql ! same, for the liquid water only (kg/m2/s)
76
77			REAL h_vcol_tot
78			! total enthalpy of vertical air column (air with water vapour,
79			! liquid and solid) (J/m2)
80
81			REAL h_dair_tot ! total enthalpy of dry air (J/m2)
82			REAL h_qw_tot ! total enthalpy of water vapour (J/m2)
83			REAL h_ql_tot ! total enthalpy of liquid water (J/m2)
84			REAL qw_tot ! total mass of water vapour (kg/m2)
85			REAL ql_tot ! total mass of liquid water (kg/m2)
86			real ec_tot ! total kinetic energy (kg/m2)
87	guez	51	REAL zairm(klon, klev) ! layer air mass (kg/m2)
88			REAL zqw_col(klon)
89			REAL zql_col(klon)
90			REAL zec_col(klon)
91			REAL zh_dair_col(klon)
92	guez	98	REAL zh_qw_col(klon), zh_ql_col(klon)
93			REAL d_h_dair, d_h_qw, d_h_ql
94	guez	51	REAL airetot, zcpvap, zcwat, zcice
95			INTEGER i, k
96			INTEGER, PARAMETER:: ndiag = 10 ! max number of diagnostic in parallel
97			integer:: pas(ndiag) = 0
98			REAL, save:: h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
99	guez	98	REAL, save:: h_ql_pre(ndiag), qw_pre(ndiag), ql_pre(ndiag)
100			REAL, save:: ec_pre(ndiag)
101	guez	51
102			!-------------------------------------------------------------
103
104			DO k = 1, klev
105			DO i = 1, klon
106			zairm(i, k) = (paprs(i, k)-paprs(i, k+1))/RG
107			ENDDO
108			END DO
109
110			! Reset variables
111			DO i = 1, klon
112			zqw_col(i)=0.
113			zql_col(i)=0.
114			zec_col(i) = 0.
115			zh_dair_col(i) = 0.
116			zh_qw_col(i) = 0.
117			zh_ql_col(i) = 0.
118			ENDDO
119
120			zcpvap=RCPV
121			zcwat=RCW
122			zcice=RCS
123
124			! Compute vertical sum for each atmospheric column
125			DO k = 1, klev
126			DO i = 1, klon
127			! Water mass
128			zqw_col(i) = zqw_col(i) + q(i, k)*zairm(i, k)
129			zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
130			! Kinetic Energy
131			zec_col(i) = zec_col(i) +0.5(u(i, k)2+v(i, k)2)zairm(i, k)
132			! Air enthalpy
133			zh_dair_col(i) = zh_dair_col(i) &
134	guez	98	+ RCPD(1.-q(i, k)-ql(i, k))zairm(i, k)*t(i, k)
135	guez	51	zh_qw_col(i) = zh_qw_col(i) + zcpvapq(i, k)zairm(i, k)*t(i, k)
136			zh_ql_col(i) = zh_ql_col(i) &
137			+ zcwatql(i, k)zairm(i, k)*t(i, k) &
138			- RLVTTql(i, k)zairm(i, k)
139			END DO
140			ENDDO
141
142			! Mean over the planet surface
143			qw_tot = 0.
144			ql_tot = 0.
145			ec_tot = 0.
146			h_vcol_tot = 0.
147			h_dair_tot = 0.
148			h_qw_tot = 0.
149			h_ql_tot = 0.
150			airetot=0.
151
152			do i=1, klon
153			qw_tot = qw_tot + zqw_col(i)*airephy(i)
154			ql_tot = ql_tot + zql_col(i)*airephy(i)
155			ec_tot = ec_tot + zec_col(i)*airephy(i)
156			h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
157			h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
158			h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
159			airetot=airetot+airephy(i)
160			END DO
161
162			qw_tot = qw_tot/airetot
163			ql_tot = ql_tot/airetot
164			ec_tot = ec_tot/airetot
165			h_dair_tot = h_dair_tot/airetot
166			h_qw_tot = h_qw_tot/airetot
167			h_ql_tot = h_ql_tot/airetot
168
169	guez	98	h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot
170	guez	51
171			! Compute the change of the atmospheric state compared to the one
172			! stored in "idiag2", and convert it in flux. This computation is
173			! performed if idiag2 /= 0 and if it is not the first call for
174			! "idiag".
175
176	guez	98	IF (idiag2 > 0 .and. pas(idiag2) /= 0) THEN
177	guez	51	d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
178			d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
179			d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
180			d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
181			d_qw = (qw_tot - qw_pre(idiag2) )/dtime
182			d_ql = (ql_tot - ql_pre(idiag2) )/dtime
183			d_ec = (ec_tot - ec_pre(idiag2) )/dtime
184	guez	98	d_qt = d_qw + d_ql
185	guez	51	ELSE
186			d_h_vcol = 0.
187			d_h_dair = 0.
188			d_h_qw = 0.
189			d_h_ql = 0.
190			d_qw = 0.
191			d_ql = 0.
192			d_ec = 0.
193			d_qt = 0.
194			ENDIF
195
196			IF (iprt >= 2) THEN
197	guez	98	print 9000, tit, pas(idiag), d_qt, d_qw, d_ql
198			print 9001, tit, pas(idiag), d_h_vcol
199			print 9002, tit, pas(idiag), d_ec
200	guez	51	END IF
201
202			! Store the new atmospheric state in "idiag"
203			pas(idiag)=pas(idiag)+1
204			h_vcol_pre(idiag) = h_vcol_tot
205			h_dair_pre(idiag) = h_dair_tot
206			h_qw_pre(idiag) = h_qw_tot
207			h_ql_pre(idiag) = h_ql_tot
208			qw_pre(idiag) = qw_tot
209			ql_pre(idiag) = ql_tot
210			ec_pre (idiag) = ec_tot
211
212	guez	98	9000 format('Physics water mass budget (kg/m2/s)', A15, 1i6, 10(1pE14.6))
213			9001 format('Physics enthalpy budget (W/m2) ', A15, 1i6, 10(F8.2))
214			9002 format('Physics kinetic energy budget (W/m2) ', A15, 1i6, 10(F8.2))
215
216	guez	51	END SUBROUTINE diagetpq
217
218			end module diagetpq_m