GCC Code Coverage Report


Directory: ./
File: dyn/leapfrog.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 138 243 56.8%
Branches: 111 218 50.9%

Line Branch Exec Source
1 !
2 ! $Id: leapfrog.F 3947 2021-06-21 15:15:02Z emillour $
3 !
4 c
5 c
6 37441 SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
7
8
9 cIM : pour sortir les param. du modele dans un fis. netcdf 110106
10 use IOIPSL
11 USE infotrac, ONLY: nqtot,ok_iso_verif
12 USE guide_mod, ONLY : guide_main
13 USE write_field, ONLY: writefield
14 USE control_mod, ONLY: nday, day_step, planet_type, offline,
15 & iconser, iphysiq, iperiod, dissip_period,
16 & iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
17 & periodav, ok_dyn_ave, output_grads_dyn
18 use exner_hyb_m, only: exner_hyb
19 use exner_milieu_m, only: exner_milieu
20 USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs
21 USE comconst_mod, ONLY: cpp, dtphys, dtvr, pi, ihf
22 USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
23 & statcl,conser,apdiss,purmats,ok_strato
24 USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
25 & start_time,dt
26
27 IMPLICIT NONE
28
29 c ...... Version du 10/01/98 ..........
30
31 c avec coordonnees verticales hybrides
32 c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
33
34 c=======================================================================
35 c
36 c Auteur: P. Le Van /L. Fairhead/F.Hourdin
37 c -------
38 c
39 c Objet:
40 c ------
41 c
42 c GCM LMD nouvelle grille
43 c
44 c=======================================================================
45 c
46 c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
47 c et possibilite d'appeler une fonction f(y) a derivee tangente
48 c hyperbolique a la place de la fonction a derivee sinusoidale.
49
50 c ... Possibilite de choisir le shema pour l'advection de
51 c q , en modifiant iadv dans traceur.def (10/02) .
52 c
53 c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
54 c Pour Van-Leer iadv=10
55 c
56 c-----------------------------------------------------------------------
57 c Declarations:
58 c -------------
59
60 include "dimensions.h"
61 include "paramet.h"
62 include "comdissnew.h"
63 include "comgeom.h"
64 include "description.h"
65 include "iniprint.h"
66 include "academic.h"
67
68 REAL,INTENT(IN) :: time_0 ! not used
69
70 c dynamical variables:
71 REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
72 REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
73 REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
74 REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
75 REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
76 REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
77 REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
78
79 REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
80 REAL pks(ip1jmp1) ! exner at the surface
81 REAL pk(ip1jmp1,llm) ! exner at mid-layer
82 REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
83 REAL phi(ip1jmp1,llm) ! geopotential
84 REAL w(ip1jmp1,llm) ! vertical velocity
85
86 real zqmin,zqmax
87
88 c variables dynamiques intermediaire pour le transport
89 REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
90
91 c variables dynamiques au pas -1
92 REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
93 REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
94 REAL massem1(ip1jmp1,llm)
95
96 c tendances dynamiques
97 REAL dv(ip1jm,llm),du(ip1jmp1,llm)
98 1 REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
99
100 c tendances de la dissipation
101 REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
102 REAL dtetadis(ip1jmp1,llm)
103
104 c tendances physiques
105 REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
106 1 REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
107
108 c variables pour le fichier histoire
109 REAL dtav ! intervalle de temps elementaire
110
111 REAL tppn(iim),tpps(iim),tpn,tps
112 c
113 INTEGER itau,itaufinp1,iav
114 ! INTEGER iday ! jour julien
115 REAL time
116
117 REAL SSUM
118 ! REAL finvmaold(ip1jmp1,llm)
119
120 cym LOGICAL lafin
121 LOGICAL :: lafin=.false.
122 INTEGER ij,iq,l
123 INTEGER ik
124
125 real time_step, t_wrt, t_ops
126
127 ! REAL rdayvrai,rdaym_ini
128 ! jD_cur: jour julien courant
129 ! jH_cur: heure julienne courante
130 REAL :: jD_cur, jH_cur
131 INTEGER :: an, mois, jour
132 REAL :: secondes
133
134 LOGICAL first,callinigrads
135 cIM : pour sortir les param. du modele dans un fis. netcdf 110106
136 save first
137 data first/.true./
138 real dt_cum
139 character*10 infile
140 integer zan, tau0, thoriid
141 integer nid_ctesGCM
142 save nid_ctesGCM
143 real degres
144 real rlong(iip1), rlatg(jjp1)
145 real zx_tmp_2d(iip1,jjp1)
146 integer ndex2d(iip1*jjp1)
147 logical ok_sync
148 parameter (ok_sync = .true.)
149 logical physic
150
151 data callinigrads/.true./
152 character*10 string10
153
154 REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
155
156 c+jld variables test conservation energie
157 REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
158 C Tendance de la temp. potentiel d (theta)/ d t due a la
159 C tansformation d'energie cinetique en energie thermique
160 C cree par la dissipation
161 REAL dtetaecdt(ip1jmp1,llm)
162 REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
163 REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
164 REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
165 CHARACTER*15 ztit
166 !IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
167 !IM SAVE ip_ebil_dyn
168 !IM DATA ip_ebil_dyn/0/
169 c-jld
170
171 character*80 dynhist_file, dynhistave_file
172 character(len=*),parameter :: modname="leapfrog"
173 character*80 abort_message
174
175 logical dissip_conservative
176 save dissip_conservative
177 data dissip_conservative/.true./
178
179 LOGICAL prem
180 save prem
181 DATA prem/.true./
182 INTEGER testita
183 PARAMETER (testita = 9)
184
185 logical , parameter :: flag_verif = .false.
186
187
188 integer itau_w ! pas de temps ecriture = itap + itau_phy
189
190
191
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 if (nday>=0) then
192 1 itaufin = nday*day_step
193 else
194 itaufin = -nday
195 endif
196 1 itaufinp1 = itaufin +1
197 1 itau = 0
198 physic=.true.
199
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if (iflag_phys==0.or.iflag_phys==2) physic=.false.
200
201 c iday = day_ini+itau/day_step
202 c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
203 c IF(time.GT.1.) THEN
204 c time = time-1.
205 c iday = iday+1
206 c ENDIF
207
208
209 c-----------------------------------------------------------------------
210 c On initialise la pression et la fonction d'Exner :
211 c --------------------------------------------------
212
213
6/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 195 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 212355 times.
✓ Branch 5 taken 195 times.
212556 dq(:,:,:)=0.
214 1 CALL pression ( ip1jmp1, ap, bp, ps, p )
215
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 if (pressure_exner) then
216 1 CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
217 else
218 CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
219 endif
220
221 c-----------------------------------------------------------------------
222 c Debut de l'integration temporelle:
223 c ----------------------------------
224
225 1 CONTINUE ! Matsuno Forward step begins here
226
227 c date: (NB: date remains unchanged for Backward step)
228 c -----
229
230 jD_cur = jD_ref + day_ini - day_ref + &
231 481 & (itau+1)/day_step
232 jH_cur = jH_ref + start_time + &
233 481 & mod(itau+1,day_step)/float(day_step)
234 481 jD_cur = jD_cur + int(jH_cur)
235 481 jH_cur = jH_cur - int(jH_cur)
236
237
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 481 times.
481 if (ok_iso_verif) then
238 call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
239 endif !if (ok_iso_verif) then
240
241
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 481 times.
481 if (ok_guide) then
242 call guide_main(itau,ucov,vcov,teta,q,masse,ps)
243 endif
244
245
246 c
247 c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
248 c CALL test_period ( ucov,vcov,teta,q,p,phis )
249 c PRINT *,' ---- Test_period apres continue OK ! -----', itau
250 c ENDIF
251 c
252
253 ! Save fields obtained at previous time step as '...m1'
254 481 CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
255 481 CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
256 481 CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
257 481 CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
258 481 CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
259
260 481 forward = .TRUE.
261 481 leapf = .FALSE.
262 481 dt = dtvr
263
264 c ... P.Le Van .26/04/94 ....
265 ! Ehouarn: finvmaold is actually not used
266 ! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
267 ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
268
269
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 481 times.
481 if (ok_iso_verif) then
270 call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
271 endif !if (ok_iso_verif) then
272
273 2 CONTINUE ! Matsuno backward or leapfrog step begins here
274
275 c-----------------------------------------------------------------------
276
277 c date: (NB: only leapfrog step requires recomputing date)
278 c -----
279
280
2/2
✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 961 times.
2881 IF (leapf) THEN
281 jD_cur = jD_ref + day_ini - day_ref +
282 1920 & (itau+1)/day_step
283 jH_cur = jH_ref + start_time +
284 1920 & mod(itau+1,day_step)/float(day_step)
285 1920 jD_cur = jD_cur + int(jH_cur)
286 1920 jH_cur = jH_cur - int(jH_cur)
287 ENDIF
288
289
290 c gestion des appels de la physique et des dissipations:
291 c ------------------------------------------------------
292 c
293 c ... P.Le Van ( 6/02/95 ) ....
294
295 2881 apphys = .FALSE.
296 2881 statcl = .FALSE.
297 2881 conser = .FALSE.
298 2881 apdiss = .FALSE.
299
300
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2881 times.
2881 IF( purmats ) THEN
301 ! Purely Matsuno time stepping
302 IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
303 IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
304 s apdiss = .TRUE.
305 IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
306 s .and. physic ) apphys = .TRUE.
307 ELSE
308 ! Leapfrog/Matsuno time stepping
309
2/2
✓ Branch 0 taken 81 times.
✓ Branch 1 taken 2800 times.
2881 IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
310
3/4
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 2401 times.
✓ Branch 2 taken 480 times.
✗ Branch 3 not taken.
2881 IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
311 480 s apdiss = .TRUE.
312
3/4
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 2401 times.
✓ Branch 2 taken 480 times.
✗ Branch 3 not taken.
2881 IF( MOD(itau+1,iphysiq).EQ.0.AND.physic ) apphys=.TRUE.
313 END IF
314
315 ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
316 ! supress dissipation step
317 if (llm.eq.1) then
318 apdiss=.false.
319 endif
320
321
322
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2881 times.
2881 if (ok_iso_verif) then
323 call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
324 endif !if (ok_iso_verif) then
325
326 c-----------------------------------------------------------------------
327 c calcul des tendances dynamiques:
328 c --------------------------------
329
330 ! compute geopotential phi()
331 2881 CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
332
333 2881 time = jD_cur + jH_cur
334 CALL caldyn
335 $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
336 2881 $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
337
338
339 c-----------------------------------------------------------------------
340 c calcul des tendances advection des traceurs (dont l'humidite)
341 c -------------------------------------------------------------
342
343
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2881 times.
2881 if (ok_iso_verif) then
344 call check_isotopes_seq(q,ip1jmp1,
345 & 'leapfrog 686: avant caladvtrac')
346 endif !if (ok_iso_verif) then
347
348
4/4
✓ Branch 0 taken 2400 times.
✓ Branch 1 taken 481 times.
✓ Branch 2 taken 1920 times.
✓ Branch 3 taken 480 times.
2881 IF( forward. OR . leapf ) THEN
349 ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
350 CALL caladvtrac(q,pbaru,pbarv,
351 * p, masse, dq, teta,
352 2401 . flxw, pk)
353 !write(*,*) 'caladvtrac 346'
354
355
356
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2401 times.
2401 IF (offline) THEN
357 Cmaf stokage du flux de masse pour traceurs OFF-LINE
358
359 CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
360 . dtvr, itau)
361
362
363 ENDIF ! of IF (offline)
364 c
365 ENDIF ! of IF( forward. OR . leapf )
366
367
368 c-----------------------------------------------------------------------
369 c integrations dynamique et traceurs:
370 c ----------------------------------
371
372
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2881 times.
2881 if (ok_iso_verif) then
373 write(*,*) 'leapfrog 720'
374 call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
375 endif !if (ok_iso_verif) then
376
377 CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
378 2881 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
379 ! $ finvmaold )
380
381
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2881 times.
2881 if (ok_iso_verif) then
382 write(*,*) 'leapfrog 724'
383 call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
384 endif !if (ok_iso_verif) then
385
386 c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
387 c
388 c-----------------------------------------------------------------------
389 c calcul des tendances physiques:
390 c -------------------------------
391 c ######## P.Le Van ( Modif le 6/02/95 ) ###########
392 c
393
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2881 times.
2881 IF( purmats ) THEN
394 IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
395 ELSE
396
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2880 times.
2881 IF( itau+1. EQ. itaufin ) lafin = .TRUE.
397 ENDIF
398 c
399 c
400
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 2401 times.
2881 IF( apphys ) THEN
401 c
402 c ....... Ajout P.Le Van ( 17/04/96 ) ...........
403 c
404
405 480 CALL pression ( ip1jmp1, ap, bp, ps, p )
406
1/2
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
480 if (pressure_exner) then
407 480 CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
408 else
409 CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
410 endif
411
412 ! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
413 ! avec dyn3dmem
414 480 CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
415
416 ! rdaym_ini = itau * dtvr / daysec
417 ! rdayvrai = rdaym_ini + day_ini
418 ! jD_cur = jD_ref + day_ini - day_ref
419 ! $ + int (itau * dtvr / daysec)
420 ! jH_cur = jH_ref + &
421 ! & (itau * dtvr / daysec - int(itau * dtvr / daysec))
422 jD_cur = jD_ref + day_ini - day_ref + &
423 480 & (itau+1)/day_step
424
425
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
480 IF (planet_type .eq."generic") THEN
426 ! AS: we make jD_cur to be pday
427 jD_cur = int(day_ini + itau/day_step)
428 ENDIF
429
430 jH_cur = jH_ref + start_time + &
431 480 & mod(itau+1,day_step)/float(day_step)
432 480 jD_cur = jD_cur + int(jH_cur)
433 480 jH_cur = jH_cur - int(jH_cur)
434 ! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
435 ! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
436 ! write(lunout,*)'current date = ',an, mois, jour, secondes
437
438 c rajout debug
439 c lafin = .true.
440
441
442 c Inbterface avec les routines de phylmd (phymars ... )
443 c -----------------------------------------------------
444
445 c+jld
446
447 c Diagnostique de conservation de l'´┐Żnergie : initialisation
448
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
480 IF (ip_ebil_dyn.ge.1 ) THEN
449 ztit='bil dyn'
450 ! Ehouarn: be careful, diagedyn is Earth-specific!
451 IF (planet_type.eq."earth") THEN
452 CALL diagedyn(ztit,2,1,1,dtphys
453 & , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
454 ENDIF
455 ENDIF ! of IF (ip_ebil_dyn.ge.1 )
456 c-jld
457 cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
458 cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
459 c IF (first) THEN
460 c first=.false.
461 c#include "ini_paramLMDZ_dyn.h"
462 c ENDIF
463 c
464 c#include "write_paramLMDZ_dyn.h"
465 c
466 ! #endif of #ifdef 1
467 CALL calfis( lafin , jD_cur, jH_cur,
468 $ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
469 $ du,dv,dteta,dq,
470 480 $ flxw,dufi,dvfi,dtetafi,dqfi,dpfi )
471 c ajout des tendances physiques:
472 c ------------------------------
473 CALL addfi( dtphys, leapf, forward ,
474 $ ucov, vcov, teta , q ,ps ,
475 480 $ dufi, dvfi, dtetafi , dqfi ,dpfi )
476 ! since addfi updates ps(), also update p(), masse() and pk()
477 480 CALL pression (ip1jmp1,ap,bp,ps,p)
478 480 CALL massdair(p,masse)
479
1/2
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
480 if (pressure_exner) then
480 480 CALL exner_hyb(ip1jmp1,ps,p,pks,pk,pkf)
481 else
482 CALL exner_milieu(ip1jmp1,ps,p,pks,pk,pkf)
483 endif
484
485
1/2
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
480 IF (ok_strato) THEN
486 480 CALL top_bound( vcov,ucov,teta,masse,dtphys)
487 ENDIF
488
489 c
490 c Diagnostique de conservation de l'´┐Żnergie : difference
491
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
480 IF (ip_ebil_dyn.ge.1 ) THEN
492 ztit='bil phys'
493 IF (planet_type.eq."earth") THEN
494 CALL diagedyn(ztit,2,1,1,dtphys
495 & , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
496 ENDIF
497 ENDIF ! of IF (ip_ebil_dyn.ge.1 )
498
499 ENDIF ! of IF( apphys )
500
501
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2881 times.
2881 IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
502 ! Academic case : Simple friction and Newtonan relaxation
503 ! -------------------------------------------------------
504 DO l=1,llm
505 DO ij=1,ip1jmp1
506 teta(ij,l)=teta(ij,l)-dtvr*
507 & (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
508 ENDDO
509 ENDDO ! of DO l=1,llm
510
511 if (planet_type.eq."giant") then
512 ! add an intrinsic heat flux at the base of the atmosphere
513 teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
514 endif
515
516 call friction(ucov,vcov,dtvr)
517
518 ! Sponge layer (if any)
519 IF (ok_strato) THEN
520 ! dufi(:,:)=0.
521 ! dvfi(:,:)=0.
522 ! dtetafi(:,:)=0.
523 ! dqfi(:,:,:)=0.
524 ! dpfi(:)=0.
525 ! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
526 CALL top_bound( vcov,ucov,teta,masse,dtvr)
527 ! CALL addfi( dtvr, leapf, forward ,
528 ! $ ucov, vcov, teta , q ,ps ,
529 ! $ dufi, dvfi, dtetafi , dqfi ,dpfi )
530 ENDIF ! of IF (ok_strato)
531 ENDIF ! of IF (iflag_phys.EQ.2)
532
533
534 c-jld
535
536 2881 CALL pression ( ip1jmp1, ap, bp, ps, p )
537
1/2
✓ Branch 0 taken 2881 times.
✗ Branch 1 not taken.
2881 if (pressure_exner) then
538 2881 CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
539 else
540 CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
541 endif
542 2881 CALL massdair(p,masse)
543
544
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2881 times.
2881 if (ok_iso_verif) then
545 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
546 endif !if (ok_iso_verif) then
547
548 c-----------------------------------------------------------------------
549 c dissipation horizontale et verticale des petites echelles:
550 c ----------------------------------------------------------
551
552
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 2401 times.
2881 IF(apdiss) THEN
553
554
555 c calcul de l'energie cinetique avant dissipation
556 480 call covcont(llm,ucov,vcov,ucont,vcont)
557 480 call enercin(vcov,ucov,vcont,ucont,ecin0)
558
559 c dissipation
560 480 CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
561
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 20386080 times.
✓ Branch 3 taken 18720 times.
20405280 ucov=ucov+dudis
562
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 19768320 times.
✓ Branch 3 taken 18720 times.
19787520 vcov=vcov+dvdis
563 c teta=teta+dtetadis
564
565
566 c------------------------------------------------------------------------
567
1/2
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
480 if (dissip_conservative) then
568 C On rajoute la tendance due a la transform. Ec -> E therm. cree
569 C lors de la dissipation
570 480 call covcont(llm,ucov,vcov,ucont,vcont)
571 480 call enercin(vcov,ucov,vcont,ucont,ecin)
572
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 20386080 times.
✓ Branch 3 taken 18720 times.
20405280 dtetaecdt= (ecin0-ecin)/ pk
573 c teta=teta+dtetaecdt
574
4/4
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 20386080 times.
✓ Branch 3 taken 18720 times.
20405280 dtetadis=dtetadis+dtetaecdt
575 endif
576
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 20386080 times.
✓ Branch 3 taken 18720 times.
20405280 teta=teta+dtetadis
577 c------------------------------------------------------------------------
578
579
580 c ....... P. Le Van ( ajout le 17/04/96 ) ...........
581 c ... Calcul de la valeur moyenne, unique de h aux poles .....
582 c
583
584
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
19200 DO l = 1, llm
585
2/2
✓ Branch 0 taken 599040 times.
✓ Branch 1 taken 18720 times.
617760 DO ij = 1,iim
586 599040 tppn(ij) = aire( ij ) * teta( ij ,l)
587 617760 tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
588 ENDDO
589 18720 tpn = SSUM(iim,tppn,1)/apoln
590 18720 tps = SSUM(iim,tpps,1)/apols
591
592
2/2
✓ Branch 0 taken 617760 times.
✓ Branch 1 taken 18720 times.
636960 DO ij = 1, iip1
593 617760 teta( ij ,l) = tpn
594 636480 teta(ij+ip1jm,l) = tps
595 ENDDO
596 ENDDO
597
598 if (1 == 0) then
599 !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
600 !!! 2) should probably not be here anyway
601 !!! but are kept for those who would want to revert to previous behaviour
602 DO ij = 1,iim
603 tppn(ij) = aire( ij ) * ps ( ij )
604 tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
605 ENDDO
606 tpn = SSUM(iim,tppn,1)/apoln
607 tps = SSUM(iim,tpps,1)/apols
608
609 DO ij = 1, iip1
610 ps( ij ) = tpn
611 ps(ij+ip1jm) = tps
612 ENDDO
613 endif ! of if (1 == 0)
614
615 END IF ! of IF(apdiss)
616
617 c ajout debug
618 c IF( lafin ) then
619 c abort_message = 'Simulation finished'
620 c call abort_gcm(modname,abort_message,0)
621 c ENDIF
622
623 c ********************************************************************
624 c ********************************************************************
625 c .... fin de l'integration dynamique et physique pour le pas itau ..
626 c ********************************************************************
627 c ********************************************************************
628
629 c preparation du pas d'integration suivant ......
630
631
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2881 times.
2881 if (ok_iso_verif) then
632 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
633 endif !if (ok_iso_verif) then
634
635
1/2
✓ Branch 0 taken 2881 times.
✗ Branch 1 not taken.
2881 IF ( .NOT.purmats ) THEN
636 c ........................................................
637 c .............. schema matsuno + leapfrog ..............
638 c ........................................................
639
640
4/4
✓ Branch 0 taken 2400 times.
✓ Branch 1 taken 481 times.
✓ Branch 2 taken 1920 times.
✓ Branch 3 taken 480 times.
2881 IF(forward. OR. leapf) THEN
641 2401 itau= itau + 1
642 c iday= day_ini+itau/day_step
643 c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
644 c IF(time.GT.1.) THEN
645 c time = time-1.
646 c iday = iday+1
647 c ENDIF
648 ENDIF
649
650
651
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2880 times.
2881 IF( itau. EQ. itaufinp1 ) then
652 if (flag_verif) then
653 write(79,*) 'ucov',ucov
654 write(80,*) 'vcov',vcov
655 write(81,*) 'teta',teta
656 write(82,*) 'ps',ps
657 write(83,*) 'q',q
658 WRITE(85,*) 'q1 = ',q(:,:,1)
659 WRITE(86,*) 'q3 = ',q(:,:,3)
660 endif
661
662 1 abort_message = 'Simulation finished'
663
664 1 call abort_gcm(modname,abort_message,0)
665 ENDIF
666 c-----------------------------------------------------------------------
667 c ecriture du fichier histoire moyenne:
668 c -------------------------------------
669
670
3/4
✓ Branch 0 taken 2400 times.
✓ Branch 1 taken 480 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2400 times.
2880 IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
671 IF(itau.EQ.itaufin) THEN
672 iav=1
673 ELSE
674 iav=0
675 ENDIF
676
677 ! ! Ehouarn: re-compute geopotential for outputs
678 480 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
679
680
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
480 IF (ok_dynzon) THEN
681 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
682 & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
683 END IF
684
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
480 IF (ok_dyn_ave) THEN
685 CALL writedynav(itau,vcov,
686 & ucov,teta,pk,phi,q,masse,ps,phis)
687 ENDIF
688
689 ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
690
691
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2880 times.
2880 if (ok_iso_verif) then
692 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
693 endif !if (ok_iso_verif) then
694
695 c-----------------------------------------------------------------------
696 c ecriture de la bande histoire:
697 c ------------------------------
698
699
1/2
✓ Branch 0 taken 2880 times.
✗ Branch 1 not taken.
2880 IF( MOD(itau,iecri).EQ.0) THEN
700 ! Ehouarn: output only during LF or Backward Matsuno
701
4/4
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 480 times.
✓ Branch 3 taken 480 times.
2880 if (leapf.or.(.not.leapf.and.(.not.forward))) then
702 2400 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
703 2400 unat=0.
704
2/2
✓ Branch 0 taken 93600 times.
✓ Branch 1 taken 2400 times.
96000 do l=1,llm
705
2/2
✓ Branch 0 taken 95752800 times.
✓ Branch 1 taken 93600 times.
95846400 unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
706
2/2
✓ Branch 0 taken 98841600 times.
✓ Branch 1 taken 93600 times.
98937600 vnat(:,l)=vcov(:,l)/cv(:)
707 enddo
708
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2400 times.
2400 if (ok_dyn_ins) then
709 ! write(lunout,*) "leapfrog: call writehist, itau=",itau
710 CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
711 ! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
712 ! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
713 ! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
714 ! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
715 ! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
716 endif ! of if (ok_dyn_ins)
717 ! For some Grads outputs of fields
718
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2400 times.
2400 if (output_grads_dyn) then
719 !
720 ! $Header$
721 !
722 if (callinigrads) then
723
724 string10='dyn'
725 call inigrads(1,iip1
726 s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
727 s ,llm,presnivs,1.
728 s ,dtvr*iperiod,string10,'dyn_zon ')
729
730 callinigrads=.false.
731
732
733 endif
734
735 string10='ps'
736 CALL wrgrads(1,1,ps,string10,string10)
737
738 string10='u'
739 CALL wrgrads(1,llm,unat,string10,string10)
740 string10='v'
741 CALL wrgrads(1,llm,vnat,string10,string10)
742 string10='teta'
743 CALL wrgrads(1,llm,teta,string10,string10)
744 do iq=1,nqtot
745 string10='q'
746 write(string10(2:2),'(i1)') iq
747 CALL wrgrads(1,llm,q(:,:,iq),string10,string10)
748 enddo
749
750 endif
751 endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
752 ENDIF ! of IF(MOD(itau,iecri).EQ.0)
753
754
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2879 times.
2880 IF(itau.EQ.itaufin) THEN
755
756
757 ! if (planet_type.eq."earth") then
758 ! Write an Earth-format restart file
759 CALL dynredem1("restart.nc",start_time,
760 1 & vcov,ucov,teta,q,masse,ps)
761 ! endif ! of if (planet_type.eq."earth")
762
763 1 CLOSE(99)
764
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if (ok_guide) then
765 ! set ok_guide to false to avoid extra output
766 ! in following forward step
767 ok_guide=.false.
768 endif
769 !!! Ehouarn: Why not stop here and now?
770 ENDIF ! of IF (itau.EQ.itaufin)
771
772 c-----------------------------------------------------------------------
773 c gestion de l'integration temporelle:
774 c ------------------------------------
775
776
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 2400 times.
2880 IF( MOD(itau,iperiod).EQ.0 ) THEN
777 GO TO 1
778
2/2
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 1440 times.
2400 ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
779
780
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 480 times.
960 IF( forward ) THEN
781 c fin du pas forward et debut du pas backward
782
783 480 forward = .FALSE.
784 480 leapf = .FALSE.
785 480 GO TO 2
786
787 ELSE
788 c fin du pas backward et debut du premier pas leapfrog
789
790 480 leapf = .TRUE.
791 480 dt = 2.*dtvr
792 480 GO TO 2
793 END IF ! of IF (forward)
794 ELSE
795
796 c ...... pas leapfrog .....
797
798 1440 leapf = .TRUE.
799 1440 dt = 2.*dtvr
800 1440 GO TO 2
801 END IF ! of IF (MOD(itau,iperiod).EQ.0)
802 ! ELSEIF (MOD(itau-1,iperiod).EQ.0)
803
804 ELSE ! of IF (.not.purmats)
805
806 if (ok_iso_verif) then
807 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
808 endif !if (ok_iso_verif) then
809
810 c ........................................................
811 c .............. schema matsuno ...............
812 c ........................................................
813 IF( forward ) THEN
814
815 itau = itau + 1
816 c iday = day_ini+itau/day_step
817 c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
818 c
819 c IF(time.GT.1.) THEN
820 c time = time-1.
821 c iday = iday+1
822 c ENDIF
823
824 forward = .FALSE.
825 IF( itau. EQ. itaufinp1 ) then
826 abort_message = 'Simulation finished'
827 call abort_gcm(modname,abort_message,0)
828 ENDIF
829 GO TO 2
830
831 ELSE ! of IF(forward) i.e. backward step
832
833 if (ok_iso_verif) then
834 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
835 endif !if (ok_iso_verif) then
836
837 IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
838 IF(itau.EQ.itaufin) THEN
839 iav=1
840 ELSE
841 iav=0
842 ENDIF
843
844 ! ! Ehouarn: re-compute geopotential for outputs
845 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
846
847 IF (ok_dynzon) THEN
848 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
849 & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
850 ENDIF
851 IF (ok_dyn_ave) THEN
852 CALL writedynav(itau,vcov,
853 & ucov,teta,pk,phi,q,masse,ps,phis)
854 ENDIF
855
856 ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
857
858 IF(MOD(itau,iecri ).EQ.0) THEN
859 c IF(MOD(itau,iecri*day_step).EQ.0) THEN
860 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
861 unat=0.
862 do l=1,llm
863 unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
864 vnat(:,l)=vcov(:,l)/cv(:)
865 enddo
866 if (ok_dyn_ins) then
867 ! write(lunout,*) "leapfrog: call writehist (b)",
868 ! & itau,iecri
869 CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
870 endif ! of if (ok_dyn_ins)
871 ! For some Grads outputs
872 if (output_grads_dyn) then
873 !
874 ! $Header$
875 !
876 if (callinigrads) then
877
878 string10='dyn'
879 call inigrads(1,iip1
880 s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
881 s ,llm,presnivs,1.
882 s ,dtvr*iperiod,string10,'dyn_zon ')
883
884 callinigrads=.false.
885
886
887 endif
888
889 string10='ps'
890 CALL wrgrads(1,1,ps,string10,string10)
891
892 string10='u'
893 CALL wrgrads(1,llm,unat,string10,string10)
894 string10='v'
895 CALL wrgrads(1,llm,vnat,string10,string10)
896 string10='teta'
897 CALL wrgrads(1,llm,teta,string10,string10)
898 do iq=1,nqtot
899 string10='q'
900 write(string10(2:2),'(i1)') iq
901 CALL wrgrads(1,llm,q(:,:,iq),string10,string10)
902 enddo
903
904 endif
905
906 ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
907
908 IF(itau.EQ.itaufin) THEN
909 ! if (planet_type.eq."earth") then
910 CALL dynredem1("restart.nc",start_time,
911 & vcov,ucov,teta,q,masse,ps)
912 ! endif ! of if (planet_type.eq."earth")
913 if (ok_guide) then
914 ! set ok_guide to false to avoid extra output
915 ! in following forward step
916 ok_guide=.false.
917 endif
918 ENDIF ! of IF(itau.EQ.itaufin)
919
920 forward = .TRUE.
921 GO TO 1
922
923 ENDIF ! of IF (forward)
924
925 END IF ! of IF(.not.purmats)
926
927 END
928