Context Navigation

leapfrog.F @ 2057

Last change on this file since 2057 was 2039, checked in by fhourdin, 11 years ago

Inclusion d'un appel supplementaire a geopot dans leapfrog dans

les versions dyn3d et dyn3dpar pour garantir la convergence
des 3 versions dyn3d/dyn3dpar/dyn3dmem.
La convergence fonctionne avec un compilateur gfortran 4.6.3 / openmpi
distribué sous ubuntu avec la nouvelle physique (NPv3.2) et le guidage.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Une seule exception : si on guide avec dyn3dpar et openMP.

Correction du guidage dans dyn3dmem.

un print supprime dans cv3_toutines.F90

Rm : On ne devrait perdre la convergence numérique avec les précédentes
versions que pour la nouvelle physique puisque le géopotentiel
n'est utilisé dans la physique que par les theermiques et sisvat.

=====================================================================

Added call to geopot in leapfrog in dyn3d and dyn3dpar in order to

insure numerical convergence with dyn3dmem.
The convergence dyn3d/dyn3dpar/dyn3dmem is satistied with
gfortran 4.6.3 / openmpi, new physics (NPv3.2) and nudging.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Only one exception : with nudging & dynd3dpar & openMP.

Bug fixing for nudging in dyn3dmem.

print supressed in cv3_routines.F90

Rm : Numerical convergence with previous releases should be lost for new
physics only since the geopotential is used only by thermals and
sisvat in the physics.

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: 26.4 KB

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

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

Download in other formats:

Original Format