Context Navigation

leapfrog.F @ 2600

Last change on this file since 2600 was 2600, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90 EM
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: 28.4 KB

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

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

Download in other formats:

Original Format