GCC Code Coverage Report


Directory: ./
File: dyn3d_common/diagedyn.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 9 0.0%
Branches: 0 2 0.0%

Line Branch Exec Source
1 !
2 ! $Id: diagedyn.F 2239 2015-03-23 07:27:30Z emillour $
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 !-----------------------------------------------------------------------
60 ! INCLUDE 'dimensions.h'
61 !
62 ! dimensions.h contient les dimensions du modele
63 ! ndm est tel que iim=2**ndm
64 !-----------------------------------------------------------------------
65
66 INTEGER iim,jjm,llm,ndm
67
68 PARAMETER (iim= 32,jjm=32,llm=39,ndm=1)
69
70 !-----------------------------------------------------------------------
71 !
72 ! $Header$
73 !
74 !
75 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
76 ! veillez n'utiliser que des ! pour les commentaires
77 ! et bien positionner les & des lignes de continuation
78 ! (les placer en colonne 6 et en colonne 73)
79 !
80 !
81 !-----------------------------------------------------------------------
82 ! INCLUDE 'paramet.h'
83
84 INTEGER iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
85 INTEGER kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
86 INTEGER ijmllm,mvar
87 INTEGER jcfil,jcfllm
88
89 PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3 &
90 & ,jjp1=jjm+1-1/jjm)
91 PARAMETER( llmp1 = llm+1, llmp2 = llm+2, llmm1 = llm-1 )
92 PARAMETER( kftd = iim/2 -ndm )
93 PARAMETER( ip1jm = iip1*jjm, ip1jmp1= iip1*jjp1 )
94 PARAMETER( ip1jmi1= ip1jm - iip1 )
95 PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
96 PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
97 PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
98
99 !-----------------------------------------------------------------------
100 !
101 ! $Header$
102 !
103 !CDK comgeom
104 COMMON/comgeom/ &
105 & cu(ip1jmp1),cv(ip1jm),unscu2(ip1jmp1),unscv2(ip1jm), &
106 & aire(ip1jmp1),airesurg(ip1jmp1),aireu(ip1jmp1), &
107 & airev(ip1jm),unsaire(ip1jmp1),apoln,apols, &
108 & unsairez(ip1jm),airuscv2(ip1jm),airvscu2(ip1jm), &
109 & aireij1(ip1jmp1),aireij2(ip1jmp1),aireij3(ip1jmp1), &
110 & aireij4(ip1jmp1),alpha1(ip1jmp1),alpha2(ip1jmp1), &
111 & alpha3(ip1jmp1),alpha4(ip1jmp1),alpha1p2(ip1jmp1), &
112 & alpha1p4(ip1jmp1),alpha2p3(ip1jmp1),alpha3p4(ip1jmp1), &
113 & fext(ip1jm),constang(ip1jmp1),rlatu(jjp1),rlatv(jjm), &
114 & rlonu(iip1),rlonv(iip1),cuvsurcv(ip1jm),cvsurcuv(ip1jm), &
115 & cvusurcu(ip1jmp1),cusurcvu(ip1jmp1),cuvscvgam1(ip1jm), &
116 & cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1), &
117 & cvuscugam2(ip1jmp1),cvscuvgam(ip1jm),cuscvugam(ip1jmp1), &
118 & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2, &
119 & unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1),unsairz_gam(ip1jm), &
120 & aivscu2gam(ip1jm),aiuscv2gam(ip1jm),xprimu(iip1),xprimv(iip1)
121
122 !
123 REAL &
124 & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,unsaire,apoln ,&
125 & apols,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4,&
126 & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 ,&
127 & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2 ,&
128 & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1,unsapolnga2&
129 & ,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2,unsairz_gam ,&
130 & aivscu2gam ,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu,cusurcvu,xprimu&
131 & , xprimv
132 !
133 !
134 ! $Header$
135 !
136 !
137 ! gestion des impressions de sorties et de d�bogage
138 ! lunout: unit� du fichier dans lequel se font les sorties
139 ! (par defaut 6, la sortie standard)
140 ! prt_level: niveau d'impression souhait� (0 = minimum)
141 !
142 INTEGER lunout, prt_level
143 COMMON /comprint/ lunout, prt_level
144
145 !#ifdef 1
146 !#include "../phylmd/YOMCST.h"
147 !#include "../phylmd/YOETHF.h"
148 !#endif
149 ! Ehouarn: for now set these parameters to what is in Earth physics...
150 ! (cf ../phylmd/suphel.h)
151 ! this should be generalized...
152 REAL,PARAMETER :: RCPD=
153 & 3.5*(1000.*(6.0221367E+23*1.380658E-23)/28.9644)
154 REAL,PARAMETER :: RCPV=
155 & 4.*(1000.*(6.0221367E+23*1.380658E-23)/18.0153)
156 REAL,PARAMETER :: RCS=RCPV
157 REAL,PARAMETER :: RCW=RCPV
158 REAL,PARAMETER :: RLSTT=2.8345E+6
159 REAL,PARAMETER :: RLVTT=2.5008E+6
160 !
161 C
162 INTEGER imjmp1
163 PARAMETER( imjmp1=iim*jjp1)
164 c Input variables
165 CHARACTER*15 tit
166 INTEGER iprt,idiag, idiag2
167 REAL dtime
168 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
169 REAL ps(ip1jmp1) ! pression au sol
170 REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
171 REAL pk (ip1jmp1,llm ) ! = (p/Pref)**kappa
172 REAL teta(ip1jmp1,llm) ! temperature potentielle
173 REAL q(ip1jmp1,llm) ! champs eau vapeur
174 REAL ql(ip1jmp1,llm) ! champs eau liquide
175
176
177 c Output variables
178 REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
179 C
180 C Local variables
181 c
182 REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
183 . , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
184 c h_vcol_tot-- total enthalpy of vertical air column
185 C (air with watter vapour, liquid and solid) (J/m2)
186 c h_dair_tot-- total enthalpy of dry air (J/m2)
187 c h_qw_tot---- total enthalpy of watter vapour (J/m2)
188 c h_ql_tot---- total enthalpy of liquid watter (J/m2)
189 c h_qs_tot---- total enthalpy of solid watter (J/m2)
190 c qw_tot------ total mass of watter vapour (kg/m2)
191 c ql_tot------ total mass of liquid watter (kg/m2)
192 c qs_tot------ total mass of solid watter (kg/m2)
193 c ec_tot------ total cinetic energy (kg/m2)
194 C
195 REAL masse(ip1jmp1,llm) ! masse d'air
196 REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
197 REAL ecin(ip1jmp1,llm)
198
199 REAL zaire(imjmp1)
200 REAL zps(imjmp1)
201 REAL zairm(imjmp1,llm)
202 REAL zecin(imjmp1,llm)
203 REAL zpaprs(imjmp1,llm)
204 REAL zpk(imjmp1,llm)
205 REAL zt(imjmp1,llm)
206 REAL zh(imjmp1,llm)
207 REAL zqw(imjmp1,llm)
208 REAL zql(imjmp1,llm)
209 REAL zqs(imjmp1,llm)
210
211 REAL zqw_col(imjmp1)
212 REAL zql_col(imjmp1)
213 REAL zqs_col(imjmp1)
214 REAL zec_col(imjmp1)
215 REAL zh_dair_col(imjmp1)
216 REAL zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1)
217 C
218 REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
219 C
220 REAL airetot, zcpvap, zcwat, zcice
221 C
222 INTEGER i, k, jj, ij , l ,ip1jjm1
223 C
224 INTEGER ndiag ! max number of diagnostic in parallel
225 PARAMETER (ndiag=10)
226 integer pas(ndiag)
227 save pas
228 data pas/ndiag*0/
229 C
230 REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
231 $ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
232 $ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
233 SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
234 $ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
235
236
237 !#ifdef 1
238 IF (planet_type=="earth") THEN
239
240 c======================================================================
241 C Compute Kinetic enrgy
242 CALL covcont ( llm , ucov , vcov , ucont, vcont )
243 CALL enercin ( vcov , ucov , vcont , ucont , ecin )
244 CALL massdair( p, masse )
245 c======================================================================
246 C
247 C
248 print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?'
249 return
250 C On ne garde les donnees que dans les colonnes i=1,iim
251 DO jj = 1,jjp1
252 ip1jjm1=iip1*(jj-1)
253 DO ij = 1,iim
254 i=iim*(jj-1)+ij
255 zaire(i)=aire(ij+ip1jjm1)
256 zps(i)=ps(ij+ip1jjm1)
257 ENDDO
258 ENDDO
259 C 3D arrays
260 DO l = 1, llm
261 DO jj = 1,jjp1
262 ip1jjm1=iip1*(jj-1)
263 DO ij = 1,iim
264 i=iim*(jj-1)+ij
265 zairm(i,l) = masse(ij+ip1jjm1,l)
266 zecin(i,l) = ecin(ij+ip1jjm1,l)
267 zpaprs(i,l) = p(ij+ip1jjm1,l)
268 zpk(i,l) = pk(ij+ip1jjm1,l)
269 zh(i,l) = teta(ij+ip1jjm1,l)
270 zqw(i,l) = q(ij+ip1jjm1,l)
271 zql(i,l) = ql(ij+ip1jjm1,l)
272 zqs(i,l) = 0.
273 ENDDO
274 ENDDO
275 ENDDO
276 C
277 C Reset variables
278 DO i = 1, imjmp1
279 zqw_col(i)=0.
280 zql_col(i)=0.
281 zqs_col(i)=0.
282 zec_col(i) = 0.
283 zh_dair_col(i) = 0.
284 zh_qw_col(i) = 0.
285 zh_ql_col(i) = 0.
286 zh_qs_col(i) = 0.
287 ENDDO
288 C
289 zcpvap=RCPV
290 zcwat=RCW
291 zcice=RCS
292 C
293 C Compute vertical sum for each atmospheric column
294 C ================================================
295 DO k = 1, llm
296 DO i = 1, imjmp1
297 C Watter mass
298 zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k)
299 zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k)
300 zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k)
301 C Cinetic Energy
302 zec_col(i) = zec_col(i)
303 $ +zecin(i,k)*zairm(i,k)
304 C Air enthalpy
305 zt(i,k)= zh(i,k) * zpk(i,k) / RCPD
306 zh_dair_col(i) = zh_dair_col(i)
307 $ + RCPD*(1.-zqw(i,k)-zql(i,k)-zqs(i,k))*zairm(i,k)*zt(i,k)
308 zh_qw_col(i) = zh_qw_col(i)
309 $ + zcpvap*zqw(i,k)*zairm(i,k)*zt(i,k)
310 zh_ql_col(i) = zh_ql_col(i)
311 $ + zcwat*zql(i,k)*zairm(i,k)*zt(i,k)
312 $ - RLVTT*zql(i,k)*zairm(i,k)
313 zh_qs_col(i) = zh_qs_col(i)
314 $ + zcice*zqs(i,k)*zairm(i,k)*zt(i,k)
315 $ - RLSTT*zqs(i,k)*zairm(i,k)
316
317 END DO
318 ENDDO
319 C
320 C Mean over the planete surface
321 C =============================
322 qw_tot = 0.
323 ql_tot = 0.
324 qs_tot = 0.
325 ec_tot = 0.
326 h_vcol_tot = 0.
327 h_dair_tot = 0.
328 h_qw_tot = 0.
329 h_ql_tot = 0.
330 h_qs_tot = 0.
331 airetot=0.
332 C
333 do i=1,imjmp1
334 qw_tot = qw_tot + zqw_col(i)
335 ql_tot = ql_tot + zql_col(i)
336 qs_tot = qs_tot + zqs_col(i)
337 ec_tot = ec_tot + zec_col(i)
338 h_dair_tot = h_dair_tot + zh_dair_col(i)
339 h_qw_tot = h_qw_tot + zh_qw_col(i)
340 h_ql_tot = h_ql_tot + zh_ql_col(i)
341 h_qs_tot = h_qs_tot + zh_qs_col(i)
342 airetot=airetot+zaire(i)
343 END DO
344 C
345 qw_tot = qw_tot/airetot
346 ql_tot = ql_tot/airetot
347 qs_tot = qs_tot/airetot
348 ec_tot = ec_tot/airetot
349 h_dair_tot = h_dair_tot/airetot
350 h_qw_tot = h_qw_tot/airetot
351 h_ql_tot = h_ql_tot/airetot
352 h_qs_tot = h_qs_tot/airetot
353 C
354 h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
355 C
356 C Compute the change of the atmospheric state compare to the one
357 C stored in "idiag2", and convert it in flux. THis computation
358 C is performed IF idiag2 /= 0 and IF it is not the first CALL
359 c for "idiag"
360 C ===================================
361 C
362 IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
363 d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
364 d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
365 d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
366 d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
367 d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
368 d_qw = (qw_tot - qw_pre(idiag2) )/dtime
369 d_ql = (ql_tot - ql_pre(idiag2) )/dtime
370 d_qs = (qs_tot - qs_pre(idiag2) )/dtime
371 d_ec = (ec_tot - ec_pre(idiag2) )/dtime
372 d_qt = d_qw + d_ql + d_qs
373 ELSE
374 d_h_vcol = 0.
375 d_h_dair = 0.
376 d_h_qw = 0.
377 d_h_ql = 0.
378 d_h_qs = 0.
379 d_qw = 0.
380 d_ql = 0.
381 d_qs = 0.
382 d_ec = 0.
383 d_qt = 0.
384 ENDIF
385 C
386 IF (iprt.ge.2) THEN
387 WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
388 9000 format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15
389 $ ,1i6,10(1pE14.6))
390 WRITE(6,9001) tit,pas(idiag), d_h_vcol
391 9001 format('Dyn3d. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
392 WRITE(6,9002) tit,pas(idiag), d_ec
393 9002 format('Dyn3d. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
394 C WRITE(6,9003) tit,pas(idiag), ec_tot
395 9003 format('Dyn3d. Cinetic Energy (W/m2) ',A15,1i6,10(E15.6))
396 WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec
397 9004 format('Dyn3d. Total Energy Budget (W/m2) ',A15,1i6,10(F8.2))
398 END IF
399 C
400 C Store the new atmospheric state in "idiag"
401 C
402 pas(idiag)=pas(idiag)+1
403 h_vcol_pre(idiag) = h_vcol_tot
404 h_dair_pre(idiag) = h_dair_tot
405 h_qw_pre(idiag) = h_qw_tot
406 h_ql_pre(idiag) = h_ql_tot
407 h_qs_pre(idiag) = h_qs_tot
408 qw_pre(idiag) = qw_tot
409 ql_pre(idiag) = ql_tot
410 qs_pre(idiag) = qs_tot
411 ec_pre (idiag) = ec_tot
412 C
413 !#else
414 ELSE
415 write(lunout,*)'diagedyn: set to function with Earth parameters'
416 ENDIF ! of if (planet_type=="earth")
417 !#endif
418 ! #endif of #ifdef 1
419 RETURN
420 END
421