Context Navigation

diagedyn.F @ 5103

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

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

Download in other formats:

Original Format