Context Navigation

leapfrog.F @ 1540

Last change on this file since 1540 was 1508, checked in by emillour, 9 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

in dynlonlat_phylonlat:
add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)

in misc:
follow up updates on wxios.F (add missing_val module variable)

in dyn3d_common:
pression.F => pression.F90
misc_mod.F90: moved from misc to dyn3d_common
added new iso_verif_dyn.F
covcont.F => covcont.F90
infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
dynetat0.F => dynetat0.F90 with some code factorization
dynredem.F => dynredem.F90 with some code factorization
added dynredem_mod.F90: routines used by dynredem
iniacademic.F90 : added isotopes-related initialization for Earth case

in dyn3d:
added check_isotopes.F
modified (isotopes) advtrac.F90, caladvtrac.F
guide_mod.F90: ported updates
leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)

in dyn3dpar:
leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
dynredem_p.F => dynredem_p.F90 and some code factorization
and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

File size: 32.1 KB

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

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

Download in other formats:

Original Format