Context Navigation

leapfrog.F @ 2276

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

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

Download in other formats:

Original Format