Context Navigation

diagedyn.F @ 5029

Last change on this file since 5029 was 4593, checked in by yann meurdesoif, 17 months ago

Replace #include (c preprocessor) by INCLUDE (fortran keyword)

in phylmd (except rrtm and ecrad) filtrez, dy3dmem and dyn3dcommon

Other directories will follow
YM

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.3 KB

Line
1	!
2	! $Id: diagedyn.F 4593 2023-06-29 13:55:54Z fhourdin $
3	!
4
5	C======================================================================
6	SUBROUTINE diagedyn(tit,iprt,idiag,idiag2,dtime
7	e , ucov , vcov , ps, p ,pk , teta , q, ql)
8	C======================================================================
9	C
10	C Purpose:
11	C Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
12	C et calcul le flux de chaleur et le flux d'eau necessaire a ces
13	C changements. Ces valeurs sont moyennees sur la surface de tout
14	C le globe et sont exprime en W/2 et kg/s/m2
15	C Outil pour diagnostiquer la conservation de l'energie
16	C et de la masse dans la dynamique.
17	C
18	C
19	c======================================================================
20	C Arguments:
21	C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
22	C iprt----input-I- PRINT level ( <=1 : no PRINT)
23	C idiag---input-I- indice dans lequel sera range les nouveaux
24	C bilans d' entalpie et de masse
25	C idiag2--input-I-les nouveaux bilans d'entalpie et de masse
26	C sont compare au bilan de d'enthalpie de masse de
27	C l'indice numero idiag2
28	C Cas parriculier : si idiag2=0, pas de comparaison, on
29	c sort directement les bilans d'enthalpie et de masse
30	C dtime----input-R- time step (s)
31	C uconv, vconv-input-R- vents covariants (m/s)
32	C ps-------input-R- Surface pressure (Pa)
33	C p--------input-R- pressure at the interfaces
34	C pk-------input-R- pk= (p/Pref)**kappa
35	c teta-----input-R- potential temperature (K)
36	c q--------input-R- vapeur d'eau (kg/kg)
37	c ql-------input-R- liquid watter (kg/kg)
38	c aire-----input-R- mesh surafce (m2)
39	c
40	C the following total value are computed by UNIT of earth surface
41	C
42	C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
43	c change (J/m2) during one time step (dtime) for the whole
44	C atmosphere (air, watter vapour, liquid and solid)
45	C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
46	C total watter (kg/m2) change during one time step (dtime),
47	C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
48	C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
49	C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
50	C
51	C
52	C J.L. Dufresne, July 2002
53	c======================================================================
54
55	USE control_mod, ONLY : planet_type
56
57	IMPLICIT NONE
58	C
59	INCLUDE "dimensions.h"
60	INCLUDE "paramet.h"
61	INCLUDE "comgeom.h"
62	INCLUDE "iniprint.h"
63
64	!#ifdef CPP_EARTH
65	! INCLUDE "../phylmd/YOMCST.h"
66	! INCLUDE "../phylmd/YOETHF.h"
67	!#endif
68	! Ehouarn: for now set these parameters to what is in Earth physics...
69	! (cf ../phylmd/suphel.h)
70	! this should be generalized...
71	REAL,PARAMETER :: RCPD=
72	& 3.5(1000.(6.0221367E+23*1.380658E-23)/28.9644)
73	REAL,PARAMETER :: RCPV=
74	& 4.(1000.(6.0221367E+23*1.380658E-23)/18.0153)
75	REAL,PARAMETER :: RCS=RCPV
76	REAL,PARAMETER :: RCW=RCPV
77	REAL,PARAMETER :: RLSTT=2.8345E+6
78	REAL,PARAMETER :: RLVTT=2.5008E+6
79	!
80	C
81	INTEGER imjmp1
82	PARAMETER( imjmp1=iim*jjp1)
83	c Input variables
84	CHARACTER*15 tit
85	INTEGER iprt,idiag, idiag2
86	REAL dtime
87	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
88	REAL ps(ip1jmp1) ! pression au sol
89	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
90	REAL pk (ip1jmp1,llm ) ! = (p/Pref)**kappa
91	REAL teta(ip1jmp1,llm) ! temperature potentielle
92	REAL q(ip1jmp1,llm) ! champs eau vapeur
93	REAL ql(ip1jmp1,llm) ! champs eau liquide
94
95
96	c Output variables
97	REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
98	C
99	C Local variables
100	c
101	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
102	. , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
103	c h_vcol_tot-- total enthalpy of vertical air column
104	C (air with watter vapour, liquid and solid) (J/m2)
105	c h_dair_tot-- total enthalpy of dry air (J/m2)
106	c h_qw_tot---- total enthalpy of watter vapour (J/m2)
107	c h_ql_tot---- total enthalpy of liquid watter (J/m2)
108	c h_qs_tot---- total enthalpy of solid watter (J/m2)
109	c qw_tot------ total mass of watter vapour (kg/m2)
110	c ql_tot------ total mass of liquid watter (kg/m2)
111	c qs_tot------ total mass of solid watter (kg/m2)
112	c ec_tot------ total cinetic energy (kg/m2)
113	C
114	REAL masse(ip1jmp1,llm) ! masse d'air
115	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
116	REAL ecin(ip1jmp1,llm)
117
118	REAL zaire(imjmp1)
119	REAL zps(imjmp1)
120	REAL zairm(imjmp1,llm)
121	REAL zecin(imjmp1,llm)
122	REAL zpaprs(imjmp1,llm)
123	REAL zpk(imjmp1,llm)
124	REAL zt(imjmp1,llm)
125	REAL zh(imjmp1,llm)
126	REAL zqw(imjmp1,llm)
127	REAL zql(imjmp1,llm)
128	REAL zqs(imjmp1,llm)
129
130	REAL zqw_col(imjmp1)
131	REAL zql_col(imjmp1)
132	REAL zqs_col(imjmp1)
133	REAL zec_col(imjmp1)
134	REAL zh_dair_col(imjmp1)
135	REAL zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1)
136	C
137	REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
138	C
139	REAL airetot, zcpvap, zcwat, zcice
140	C
141	INTEGER i, k, jj, ij , l ,ip1jjm1
142	C
143	INTEGER ndiag ! max number of diagnostic in parallel
144	PARAMETER (ndiag=10)
145	integer pas(ndiag)
146	save pas
147	data pas/ndiag*0/
148	C
149	REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
150	$ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
151	$ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
152	SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
153	$ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
154
155
156	!#ifdef CPP_EARTH
157	IF (planet_type=="earth") THEN
158
159	c======================================================================
160	C Compute Kinetic enrgy
161	CALL covcont ( llm , ucov , vcov , ucont, vcont )
162	CALL enercin ( vcov , ucov , vcont , ucont , ecin )
163	CALL massdair( p, masse )
164	c======================================================================
165	C
166	C
167	print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?'
168	return
169	C On ne garde les donnees que dans les colonnes i=1,iim
170	DO jj = 1,jjp1
171	ip1jjm1=iip1*(jj-1)
172	DO ij = 1,iim
173	i=iim*(jj-1)+ij
174	zaire(i)=aire(ij+ip1jjm1)
175	zps(i)=ps(ij+ip1jjm1)
176	ENDDO
177	ENDDO
178	C 3D arrays
179	DO l = 1, llm
180	DO jj = 1,jjp1
181	ip1jjm1=iip1*(jj-1)
182	DO ij = 1,iim
183	i=iim*(jj-1)+ij
184	zairm(i,l) = masse(ij+ip1jjm1,l)
185	zecin(i,l) = ecin(ij+ip1jjm1,l)
186	zpaprs(i,l) = p(ij+ip1jjm1,l)
187	zpk(i,l) = pk(ij+ip1jjm1,l)
188	zh(i,l) = teta(ij+ip1jjm1,l)
189	zqw(i,l) = q(ij+ip1jjm1,l)
190	zql(i,l) = ql(ij+ip1jjm1,l)
191	zqs(i,l) = 0.
192	ENDDO
193	ENDDO
194	ENDDO
195	C
196	C Reset variables
197	DO i = 1, imjmp1
198	zqw_col(i)=0.
199	zql_col(i)=0.
200	zqs_col(i)=0.
201	zec_col(i) = 0.
202	zh_dair_col(i) = 0.
203	zh_qw_col(i) = 0.
204	zh_ql_col(i) = 0.
205	zh_qs_col(i) = 0.
206	ENDDO
207	C
208	zcpvap=RCPV
209	zcwat=RCW
210	zcice=RCS
211	C
212	C Compute vertical sum for each atmospheric column
213	C ================================================
214	DO k = 1, llm
215	DO i = 1, imjmp1
216	C Watter mass
217	zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k)
218	zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k)
219	zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k)
220	C Cinetic Energy
221	zec_col(i) = zec_col(i)
222	$ +zecin(i,k)*zairm(i,k)
223	C Air enthalpy
224	zt(i,k)= zh(i,k) * zpk(i,k) / RCPD
225	zh_dair_col(i) = zh_dair_col(i)
226	$ + RCPD(1.-zqw(i,k)-zql(i,k)-zqs(i,k))zairm(i,k)*zt(i,k)
227	zh_qw_col(i) = zh_qw_col(i)
228	$ + zcpvapzqw(i,k)zairm(i,k)*zt(i,k)
229	zh_ql_col(i) = zh_ql_col(i)
230	$ + zcwatzql(i,k)zairm(i,k)*zt(i,k)
231	$ - RLVTTzql(i,k)zairm(i,k)
232	zh_qs_col(i) = zh_qs_col(i)
233	$ + zcicezqs(i,k)zairm(i,k)*zt(i,k)
234	$ - RLSTTzqs(i,k)zairm(i,k)
235
236	END DO
237	ENDDO
238	C
239	C Mean over the planete surface
240	C =============================
241	qw_tot = 0.
242	ql_tot = 0.
243	qs_tot = 0.
244	ec_tot = 0.
245	h_vcol_tot = 0.
246	h_dair_tot = 0.
247	h_qw_tot = 0.
248	h_ql_tot = 0.
249	h_qs_tot = 0.
250	airetot=0.
251	C
252	do i=1,imjmp1
253	qw_tot = qw_tot + zqw_col(i)
254	ql_tot = ql_tot + zql_col(i)
255	qs_tot = qs_tot + zqs_col(i)
256	ec_tot = ec_tot + zec_col(i)
257	h_dair_tot = h_dair_tot + zh_dair_col(i)
258	h_qw_tot = h_qw_tot + zh_qw_col(i)
259	h_ql_tot = h_ql_tot + zh_ql_col(i)
260	h_qs_tot = h_qs_tot + zh_qs_col(i)
261	airetot=airetot+zaire(i)
262	END DO
263	C
264	qw_tot = qw_tot/airetot
265	ql_tot = ql_tot/airetot
266	qs_tot = qs_tot/airetot
267	ec_tot = ec_tot/airetot
268	h_dair_tot = h_dair_tot/airetot
269	h_qw_tot = h_qw_tot/airetot
270	h_ql_tot = h_ql_tot/airetot
271	h_qs_tot = h_qs_tot/airetot
272	C
273	h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
274	C
275	C Compute the change of the atmospheric state compare to the one
276	C stored in "idiag2", and convert it in flux. THis computation
277	C is performed IF idiag2 /= 0 and IF it is not the first CALL
278	c for "idiag"
279	C ===================================
280	C
281	IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
282	d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
283	d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
284	d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
285	d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
286	d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
287	d_qw = (qw_tot - qw_pre(idiag2) )/dtime
288	d_ql = (ql_tot - ql_pre(idiag2) )/dtime
289	d_qs = (qs_tot - qs_pre(idiag2) )/dtime
290	d_ec = (ec_tot - ec_pre(idiag2) )/dtime
291	d_qt = d_qw + d_ql + d_qs
292	ELSE
293	d_h_vcol = 0.
294	d_h_dair = 0.
295	d_h_qw = 0.
296	d_h_ql = 0.
297	d_h_qs = 0.
298	d_qw = 0.
299	d_ql = 0.
300	d_qs = 0.
301	d_ec = 0.
302	d_qt = 0.
303	ENDIF
304	C
305	IF (iprt.ge.2) THEN
306	WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
307	9000 format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15
308	$ ,1i6,10(1pE14.6))
309	WRITE(6,9001) tit,pas(idiag), d_h_vcol
310	9001 format('Dyn3d. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
311	WRITE(6,9002) tit,pas(idiag), d_ec
312	9002 format('Dyn3d. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
313	C WRITE(6,9003) tit,pas(idiag), ec_tot
314	9003 format('Dyn3d. Cinetic Energy (W/m2) ',A15,1i6,10(E15.6))
315	WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec
316	9004 format('Dyn3d. Total Energy Budget (W/m2) ',A15,1i6,10(F8.2))
317	END IF
318	C
319	C Store the new atmospheric state in "idiag"
320	C
321	pas(idiag)=pas(idiag)+1
322	h_vcol_pre(idiag) = h_vcol_tot
323	h_dair_pre(idiag) = h_dair_tot
324	h_qw_pre(idiag) = h_qw_tot
325	h_ql_pre(idiag) = h_ql_tot
326	h_qs_pre(idiag) = h_qs_tot
327	qw_pre(idiag) = qw_tot
328	ql_pre(idiag) = ql_tot
329	qs_pre(idiag) = qs_tot
330	ec_pre (idiag) = ec_tot
331	C
332	!#else
333	ELSE
334	write(lunout,*)'diagedyn: set to function with Earth parameters'
335	ENDIF ! of if (planet_type=="earth")
336	!#endif
337	! #endif of #ifdef CPP_EARTH
338	RETURN
339	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format