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