Context Navigation

diagedyn.F @ 5089

Last change on this file since 5089 was 5089, checked in by abarral, 2 months ago
Remove unused (2015+) CPP_EARTH Start refactoring makelmdz_fcm to prepare investigating fcm upgrade
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.2 KB

Line
1	!
2	! $Id: diagedyn.F 5089 2024-07-20 15:25:07Z abarral $
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	! Ehouarn: for now set these parameters to what is in Earth physics...
65	! (cf ../phylmd/suphel.h)
66	! this should be generalized...
67	REAL,PARAMETER :: RCPD=
68	& 3.5(1000.(6.0221367E+23*1.380658E-23)/28.9644)
69	REAL,PARAMETER :: RCPV=
70	& 4.(1000.(6.0221367E+23*1.380658E-23)/18.0153)
71	REAL,PARAMETER :: RCS=RCPV
72	REAL,PARAMETER :: RCW=RCPV
73	REAL,PARAMETER :: RLSTT=2.8345E+6
74	REAL,PARAMETER :: RLVTT=2.5008E+6
75	!
76	C
77	INTEGER imjmp1
78	PARAMETER( imjmp1=iim*jjp1)
79	c Input variables
80	CHARACTER*15 tit
81	INTEGER iprt,idiag, idiag2
82	REAL dtime
83	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
84	REAL ps(ip1jmp1) ! pression au sol
85	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
86	REAL pk (ip1jmp1,llm ) ! = (p/Pref)**kappa
87	REAL teta(ip1jmp1,llm) ! temperature potentielle
88	REAL q(ip1jmp1,llm) ! champs eau vapeur
89	REAL ql(ip1jmp1,llm) ! champs eau liquide
90
91
92	c Output variables
93	REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
94	C
95	C Local variables
96	c
97	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
98	. , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
99	c h_vcol_tot-- total enthalpy of vertical air column
100	C (air with watter vapour, liquid and solid) (J/m2)
101	c h_dair_tot-- total enthalpy of dry air (J/m2)
102	c h_qw_tot---- total enthalpy of watter vapour (J/m2)
103	c h_ql_tot---- total enthalpy of liquid watter (J/m2)
104	c h_qs_tot---- total enthalpy of solid watter (J/m2)
105	c qw_tot------ total mass of watter vapour (kg/m2)
106	c ql_tot------ total mass of liquid watter (kg/m2)
107	c qs_tot------ total mass of solid watter (kg/m2)
108	c ec_tot------ total cinetic energy (kg/m2)
109	C
110	REAL masse(ip1jmp1,llm) ! masse d'air
111	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
112	REAL ecin(ip1jmp1,llm)
113
114	REAL zaire(imjmp1)
115	REAL zps(imjmp1)
116	REAL zairm(imjmp1,llm)
117	REAL zecin(imjmp1,llm)
118	REAL zpaprs(imjmp1,llm)
119	REAL zpk(imjmp1,llm)
120	REAL zt(imjmp1,llm)
121	REAL zh(imjmp1,llm)
122	REAL zqw(imjmp1,llm)
123	REAL zql(imjmp1,llm)
124	REAL zqs(imjmp1,llm)
125
126	REAL zqw_col(imjmp1)
127	REAL zql_col(imjmp1)
128	REAL zqs_col(imjmp1)
129	REAL zec_col(imjmp1)
130	REAL zh_dair_col(imjmp1)
131	REAL zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1)
132	C
133	REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
134	C
135	REAL airetot, zcpvap, zcwat, zcice
136	C
137	INTEGER i, k, jj, ij , l ,ip1jjm1
138	C
139	INTEGER ndiag ! max number of diagnostic in parallel
140	PARAMETER (ndiag=10)
141	integer pas(ndiag)
142	save pas
143	data pas/ndiag*0/
144	C
145	REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
146	$ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
147	$ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
148	SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
149	$ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
150
151
152	IF (planet_type=="earth") THEN
153
154	c======================================================================
155	C Compute Kinetic enrgy
156	CALL covcont ( llm , ucov , vcov , ucont, vcont )
157	CALL enercin ( vcov , ucov , vcont , ucont , ecin )
158	CALL massdair( p, masse )
159	c======================================================================
160	C
161	C
162	print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?'
163	return
164	C On ne garde les donnees que dans les colonnes i=1,iim
165	DO jj = 1,jjp1
166	ip1jjm1=iip1*(jj-1)
167	DO ij = 1,iim
168	i=iim*(jj-1)+ij
169	zaire(i)=aire(ij+ip1jjm1)
170	zps(i)=ps(ij+ip1jjm1)
171	ENDDO
172	ENDDO
173	C 3D arrays
174	DO l = 1, llm
175	DO jj = 1,jjp1
176	ip1jjm1=iip1*(jj-1)
177	DO ij = 1,iim
178	i=iim*(jj-1)+ij
179	zairm(i,l) = masse(ij+ip1jjm1,l)
180	zecin(i,l) = ecin(ij+ip1jjm1,l)
181	zpaprs(i,l) = p(ij+ip1jjm1,l)
182	zpk(i,l) = pk(ij+ip1jjm1,l)
183	zh(i,l) = teta(ij+ip1jjm1,l)
184	zqw(i,l) = q(ij+ip1jjm1,l)
185	zql(i,l) = ql(ij+ip1jjm1,l)
186	zqs(i,l) = 0.
187	ENDDO
188	ENDDO
189	ENDDO
190	C
191	C Reset variables
192	DO i = 1, imjmp1
193	zqw_col(i)=0.
194	zql_col(i)=0.
195	zqs_col(i)=0.
196	zec_col(i) = 0.
197	zh_dair_col(i) = 0.
198	zh_qw_col(i) = 0.
199	zh_ql_col(i) = 0.
200	zh_qs_col(i) = 0.
201	ENDDO
202	C
203	zcpvap=RCPV
204	zcwat=RCW
205	zcice=RCS
206	C
207	C Compute vertical sum for each atmospheric column
208	C ================================================
209	DO k = 1, llm
210	DO i = 1, imjmp1
211	C Watter mass
212	zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k)
213	zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k)
214	zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k)
215	C Cinetic Energy
216	zec_col(i) = zec_col(i)
217	$ +zecin(i,k)*zairm(i,k)
218	C Air enthalpy
219	zt(i,k)= zh(i,k) * zpk(i,k) / RCPD
220	zh_dair_col(i) = zh_dair_col(i)
221	$ + RCPD(1.-zqw(i,k)-zql(i,k)-zqs(i,k))zairm(i,k)*zt(i,k)
222	zh_qw_col(i) = zh_qw_col(i)
223	$ + zcpvapzqw(i,k)zairm(i,k)*zt(i,k)
224	zh_ql_col(i) = zh_ql_col(i)
225	$ + zcwatzql(i,k)zairm(i,k)*zt(i,k)
226	$ - RLVTTzql(i,k)zairm(i,k)
227	zh_qs_col(i) = zh_qs_col(i)
228	$ + zcicezqs(i,k)zairm(i,k)*zt(i,k)
229	$ - RLSTTzqs(i,k)zairm(i,k)
230
231	END DO
232	ENDDO
233	C
234	C Mean over the planete surface
235	C =============================
236	qw_tot = 0.
237	ql_tot = 0.
238	qs_tot = 0.
239	ec_tot = 0.
240	h_vcol_tot = 0.
241	h_dair_tot = 0.
242	h_qw_tot = 0.
243	h_ql_tot = 0.
244	h_qs_tot = 0.
245	airetot=0.
246	C
247	do i=1,imjmp1
248	qw_tot = qw_tot + zqw_col(i)
249	ql_tot = ql_tot + zql_col(i)
250	qs_tot = qs_tot + zqs_col(i)
251	ec_tot = ec_tot + zec_col(i)
252	h_dair_tot = h_dair_tot + zh_dair_col(i)
253	h_qw_tot = h_qw_tot + zh_qw_col(i)
254	h_ql_tot = h_ql_tot + zh_ql_col(i)
255	h_qs_tot = h_qs_tot + zh_qs_col(i)
256	airetot=airetot+zaire(i)
257	END DO
258	C
259	qw_tot = qw_tot/airetot
260	ql_tot = ql_tot/airetot
261	qs_tot = qs_tot/airetot
262	ec_tot = ec_tot/airetot
263	h_dair_tot = h_dair_tot/airetot
264	h_qw_tot = h_qw_tot/airetot
265	h_ql_tot = h_ql_tot/airetot
266	h_qs_tot = h_qs_tot/airetot
267	C
268	h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
269	C
270	C Compute the change of the atmospheric state compare to the one
271	C stored in "idiag2", and convert it in flux. THis computation
272	C is performed IF idiag2 /= 0 and IF it is not the first CALL
273	c for "idiag"
274	C ===================================
275	C
276	IF ( (idiag2>0) .and. (pas(idiag2) /= 0) ) THEN
277	d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
278	d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
279	d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
280	d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
281	d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
282	d_qw = (qw_tot - qw_pre(idiag2) )/dtime
283	d_ql = (ql_tot - ql_pre(idiag2) )/dtime
284	d_qs = (qs_tot - qs_pre(idiag2) )/dtime
285	d_ec = (ec_tot - ec_pre(idiag2) )/dtime
286	d_qt = d_qw + d_ql + d_qs
287	ELSE
288	d_h_vcol = 0.
289	d_h_dair = 0.
290	d_h_qw = 0.
291	d_h_ql = 0.
292	d_h_qs = 0.
293	d_qw = 0.
294	d_ql = 0.
295	d_qs = 0.
296	d_ec = 0.
297	d_qt = 0.
298	ENDIF
299	C
300	IF (iprt>=2) THEN
301	WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
302	9000 format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15
303	$ ,1i6,10(1pE14.6))
304	WRITE(6,9001) tit,pas(idiag), d_h_vcol
305	9001 format('Dyn3d. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
306	WRITE(6,9002) tit,pas(idiag), d_ec
307	9002 format('Dyn3d. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
308	C WRITE(6,9003) tit,pas(idiag), ec_tot
309	9003 format('Dyn3d. Cinetic Energy (W/m2) ',A15,1i6,10(E15.6))
310	WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec
311	9004 format('Dyn3d. Total Energy Budget (W/m2) ',A15,1i6,10(F8.2))
312	END IF
313	C
314	C Store the new atmospheric state in "idiag"
315	C
316	pas(idiag)=pas(idiag)+1
317	h_vcol_pre(idiag) = h_vcol_tot
318	h_dair_pre(idiag) = h_dair_tot
319	h_qw_pre(idiag) = h_qw_tot
320	h_ql_pre(idiag) = h_ql_tot
321	h_qs_pre(idiag) = h_qs_tot
322	qw_pre(idiag) = qw_tot
323	ql_pre(idiag) = ql_tot
324	qs_pre(idiag) = qs_tot
325	ec_pre (idiag) = ec_tot
326	C
327	!#else
328	ELSE
329	write(lunout,*)'diagedyn: set to function with Earth parameters'
330	ENDIF ! of if (planet_type=="earth")
331	!#endif
332	RETURN
333	END

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

Download in other formats:

Original Format