Context Navigation

leapfrog.F @ 5095

Last change on this file since 5095 was 5091, checked in by abarral, 4 months ago
Move lmdz_netcdf_format.F90 -> lmdz_cppkeys_wrapper.F90 to handle other CPP keys Replace all (except wrapper) use of CPP_PHYS by fortran logical Refactor makelmdz_fcm (put blocks into functions, use modern bash)
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.2 KB

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

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

Download in other formats:

Original Format