Context Navigation

leapfrog.F @ 1824

Last change on this file since 1824 was 1824, checked in by emillour, 7 years ago

Common dynamics:

enable possiblity to store multiple time steps in the restart.nc file (flag "ecritstart" gives the frequency, in dynamical steps).
fixed dynredem_mod.F90 to correctly write multiple time steps.
fixed computation of JH_cur in the mars case where "hour_ini" contains the initial time of day read from the start.nc file
minor fix in dynetat0.F90

File size: 34.6 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	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
540
541	jD_cur = jD_ref + day_ini - day_ref + &
542	& (itau+1)/day_step
543
544	IF ((planet_type .eq."generic").or.
545	& (planet_type .eq."mars")) THEN
546	! AS: we make jD_cur to be pday
547	jD_cur = int(day_ini + itau/day_step)
548	ENDIF
549
550	IF (planet_type .eq. "mars") THEN
551	jH_cur = jH_ref + hour_ini + &
552	& mod(itau,day_step)/float(day_step)
553	ELSE IF (planet_type .eq. "generic") THEN
554	jH_cur = jH_ref + start_time + &
555	& mod(itau,day_step)/float(day_step)
556	ELSE
557	jH_cur = jH_ref + start_time + &
558	& mod(itau+1,day_step)/float(day_step)
559	ENDIF
560	jD_cur = jD_cur + int(jH_cur)
561	jH_cur = jH_cur - int(jH_cur)
562	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
563	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
564	! write(lunout,*)'current date = ',an, mois, jour, secondes
565
566	c rajout debug
567	c lafin = .true.
568
569
570	c Interface avec les routines de phylmd (phymars ... )
571	c -----------------------------------------------------
572
573	c+jld
574
575	c Diagnostique de conservation de l'énergie : initialisation
576	IF (ip_ebil_dyn.ge.1 ) THEN
577	ztit='bil dyn'
578	! Ehouarn: be careful, diagedyn is Earth-specific!
579	IF (planet_type.eq."earth") THEN
580	CALL diagedyn(ztit,2,1,1,dtphys
581	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
582	ENDIF
583	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
584	c-jld
585	#ifdef CPP_IOIPSL
586	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
587	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
588	c IF (first) THEN
589	c first=.false.
590	c#include "ini_paramLMDZ_dyn.h"
591	c ENDIF
592	c
593	c#include "write_paramLMDZ_dyn.h"
594	c
595	#endif
596	! #endif of #ifdef CPP_IOIPSL
597
598	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
599
600	CALL calfis( lafin , jD_cur, jH_cur,
601	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
602	$ du,dv,dteta,dq,
603	$ flxw,
604	$ dufi,dvfi,dtetafi,dqfi,dpfi )
605
606	c call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
607	c call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
608	c call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
609
610	c ajout des tendances physiques:
611	c ------------------------------
612	CALL addfi( dtphys, leapf, forward ,
613	$ ucov, vcov, teta , q ,ps ,
614	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
615	! since addfi updates ps(), also update p(), masse() and pk()
616	CALL pression (ip1jmp1,ap,bp,ps,p)
617	CALL massdair(p,masse)
618	if (pressure_exner) then
619	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
620	else
621	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
622	endif
623
624	c Couche superieure :
625	c -------------------
626	IF (iflag_top_bound > 0) THEN
627	CALL top_bound(vcov,ucov,teta,masse,dtphys,dutop)
628	dutop(:,:)=dutop(:,:)/dtphys ! convert to a tendency in (m/s)/s
629	ENDIF
630
631	c Diagnostique de conservation de l'énergie : difference
632	IF (ip_ebil_dyn.ge.1 ) THEN
633	ztit='bil phys'
634	IF (planet_type.eq."earth") THEN
635	CALL diagedyn(ztit,2,1,1,dtphys
636	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
637	ENDIF
638	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
639
640	ENDIF ! of IF( apphys )
641
642	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
643	! Academic case : Simple friction and Newtonan relaxation
644	! -------------------------------------------------------
645	DO l=1,llm
646	DO ij=1,ip1jmp1
647	teta(ij,l)=teta(ij,l)-dtvr*
648	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
649	ENDDO
650	ENDDO ! of DO l=1,llm
651
652	if (planet_type.eq."giant") then
653	! add an intrinsic heat flux at the base of the atmosphere
654	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
655	endif
656
657	call friction(ucov,vcov,dtvr)
658
659
660	! Sponge layer (if any)
661	IF (ok_strato) THEN
662	CALL top_bound(vcov,ucov,teta,masse,dtvr,dutop)
663	dutop(:,:)=dutop(:,:)/dtvr ! convert to a tendency in (m/s)/s
664	ENDIF ! of IF (ok_strato)
665	ENDIF ! of IF (iflag_phys.EQ.2)
666
667
668	c-jld
669
670	CALL pression ( ip1jmp1, ap, bp, ps, p )
671	if (pressure_exner) then
672	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
673	else
674	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
675	endif
676	CALL massdair(p,masse)
677
678	if (ok_iso_verif) then
679	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
680	endif !if (ok_iso_verif) then
681
682	c-----------------------------------------------------------------------
683	c dissipation horizontale et verticale des petites echelles:
684	c ----------------------------------------------------------
685
686	IF(apdiss) THEN
687
688	! sponge layer
689	if (callsponge) then
690	CALL sponge(ucov,vcov,teta,ps,dtdiss,mode_sponge)
691	endif
692
693	c calcul de l'energie cinetique avant dissipation
694	call covcont(llm,ucov,vcov,ucont,vcont)
695	call enercin(vcov,ucov,vcont,ucont,ecin0)
696
697	c dissipation
698	! ADAPTATION GCM POUR CP(T)
699	call tpot2t(ijp1llm,teta,temp,pk)
700
701	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
702	ucov=ucov+dudis
703	vcov=vcov+dvdis
704	dudis=dudis/dtdiss ! passage en (m/s)/s
705	dvdis=dvdis/dtdiss ! passage en (m/s)/s
706
707	c------------------------------------------------------------------------
708	if (dissip_conservative) then
709	C On rajoute la tendance due a la transform. Ec -> E therm. cree
710	C lors de la dissipation
711	call covcont(llm,ucov,vcov,ucont,vcont)
712	call enercin(vcov,ucov,vcont,ucont,ecin)
713	! ADAPTATION GCM POUR CP(T)
714	do ij=1,ip1jmp1
715	do l=1,llm
716	dtec = (ecin0(ij,l)-ecin(ij,l))/cpdet(temp(ij,l))
717	temp(ij,l) = temp(ij,l) + dtec
718	enddo
719	enddo
720	call t2tpot(ijp1llm,temp,ztetaec,pk)
721	dtetaecdt=ztetaec-teta
722	dtetadis=dtetadis+dtetaecdt
723	endif
724	teta=teta+dtetadis
725	dtetadis=dtetadis/dtdiss ! passage en K/s
726	c------------------------------------------------------------------------
727
728
729	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
730	c ... Calcul de la valeur moyenne, unique de h aux poles .....
731	c
732
733	DO l = 1, llm
734	DO ij = 1,iim
735	tppn(ij) = aire( ij ) * teta( ij ,l)
736	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
737	ENDDO
738	tpn = SSUM(iim,tppn,1)/apoln
739	tps = SSUM(iim,tpps,1)/apols
740
741	DO ij = 1, iip1
742	teta( ij ,l) = tpn
743	teta(ij+ip1jm,l) = tps
744	ENDDO
745	ENDDO
746
747	if (1 == 0) then
748	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
749	!!! 2) should probably not be here anyway
750	!!! but are kept for those who would want to revert to previous behaviour
751	DO ij = 1,iim
752	tppn(ij) = aire( ij ) * ps ( ij )
753	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
754	ENDDO
755	tpn = SSUM(iim,tppn,1)/apoln
756	tps = SSUM(iim,tpps,1)/apols
757
758	DO ij = 1, iip1
759	ps( ij ) = tpn
760	ps(ij+ip1jm) = tps
761	ENDDO
762	endif ! of if (1 == 0)
763
764	END IF ! of IF(apdiss)
765
766	c ajout debug
767	c IF( lafin ) then
768	c abort_message = 'Simulation finished'
769	c call abort_gcm(modname,abort_message,0)
770	c ENDIF
771
772	c ********************************************************************
773	c ********************************************************************
774	c .... fin de l'integration dynamique et physique pour le pas itau ..
775	c ********************************************************************
776	c ********************************************************************
777
778	c preparation du pas d'integration suivant ......
779
780	if (ok_iso_verif) then
781	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
782	endif !if (ok_iso_verif) then
783
784	IF ( .NOT.purmats ) THEN
785	c ........................................................
786	c .............. schema matsuno + leapfrog ..............
787	c ........................................................
788
789	IF(forward. OR. leapf) THEN
790	itau= itau + 1
791	c iday= day_ini+itau/day_step
792	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
793	c IF(time.GT.1.) THEN
794	c time = time-1.
795	c iday = iday+1
796	c ENDIF
797	ENDIF
798
799
800	IF( itau. EQ. itaufinp1 ) then
801	if (flag_verif) then
802	write(79,*) 'ucov',ucov
803	write(80,*) 'vcov',vcov
804	write(81,*) 'teta',teta
805	write(82,*) 'ps',ps
806	write(83,*) 'q',q
807	WRITE(85,*) 'q1 = ',q(:,:,1)
808	WRITE(86,*) 'q3 = ',q(:,:,3)
809	endif
810
811	abort_message = 'Simulation finished'
812
813	call abort_gcm(modname,abort_message,0)
814	ENDIF
815	c-----------------------------------------------------------------------
816	c ecriture du fichier histoire moyenne:
817	c -------------------------------------
818
819	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
820	IF(itau.EQ.itaufin) THEN
821	iav=1
822	ELSE
823	iav=0
824	ENDIF
825
826	! ! Ehouarn: re-compute geopotential for outputs
827	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
828
829	IF (ok_dynzon) THEN
830	#ifdef CPP_IOIPSL
831	c les traceurs ne sont pas sortis, trop lourd.
832	c Peut changer eventuellement si besoin.
833	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
834	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
835	& du,dudis,dutop,dufi)
836	#endif
837	END IF
838	IF (ok_dyn_ave) THEN
839	#ifdef CPP_IOIPSL
840	CALL writedynav(itau,vcov,
841	& ucov,teta,pk,phi,q,masse,ps,phis)
842	#endif
843	ENDIF
844
845	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
846
847	if (ok_iso_verif) then
848	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
849	endif !if (ok_iso_verif) then
850
851	c-----------------------------------------------------------------------
852	c ecriture de la bande histoire:
853	c ------------------------------
854
855	IF( MOD(itau,iecri).EQ.0) THEN
856	! Ehouarn: output only during LF or Backward Matsuno
857	if (leapf.or.(.not.leapf.and.(.not.forward))) then
858	! ADAPTATION GCM POUR CP(T)
859	call tpot2t(ijp1llm,teta,temp,pk)
860	tsurpk = cpp*temp/pk
861	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
862	unat=0.
863	do l=1,llm
864	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
865	vnat(:,l)=vcov(:,l)/cv(:)
866	enddo
867	#ifdef CPP_IOIPSL
868	if (ok_dyn_ins) then
869	! write(lunout,*) "leapfrog: call writehist, itau=",itau
870	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
871	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
872	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
873	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
874	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
875	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
876	endif ! of if (ok_dyn_ins)
877	#endif
878	! For some Grads outputs of fields
879	if (output_grads_dyn) then
880	#include "write_grads_dyn.h"
881	endif
882	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
883	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
884
885	c Determine whether to write to the restart.nc file
886	c Decision can't be made in one IF statement as if
887	c ecritstart==0 there will be a divide-by-zero error
888	lrestart = .false.
889	IF (itau.EQ.itaufin) THEN
890	lrestart = .true.
891	ELSE IF (ecritstart.GT.0) THEN
892	IF (MOD(itau,ecritstart).EQ.0) lrestart = .true.
893	ENDIF
894
895	c Write to the restart.nc file
896	IF (lrestart) THEN
897	if (planet_type=="mars") then
898	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
899	& vcov,ucov,teta,q,masse,ps)
900	else
901	CALL dynredem1("restart.nc",start_time,
902	& vcov,ucov,teta,q,masse,ps)
903	endif
904	CLOSE(99)
905	!!! Ehouarn: Why not stop here and now?
906	ENDIF ! of IF (lrestart)
907
908	c-----------------------------------------------------------------------
909	c gestion de l'integration temporelle:
910	c ------------------------------------
911
912	IF( MOD(itau,iperiod).EQ.0 ) THEN
913	GO TO 1
914	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
915
916	IF( forward ) THEN
917	c fin du pas forward et debut du pas backward
918
919	forward = .FALSE.
920	leapf = .FALSE.
921	GO TO 2
922
923	ELSE
924	c fin du pas backward et debut du premier pas leapfrog
925
926	leapf = .TRUE.
927	dt = 2.*dtvr
928	GO TO 2
929	END IF ! of IF (forward)
930	ELSE
931
932	c ...... pas leapfrog .....
933
934	leapf = .TRUE.
935	dt = 2.*dtvr
936	GO TO 2
937	END IF ! of IF (MOD(itau,iperiod).EQ.0)
938	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
939
940	ELSE ! of IF (.not.purmats)
941
942	if (ok_iso_verif) then
943	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
944	endif !if (ok_iso_verif) then
945
946	c ........................................................
947	c .............. schema matsuno ...............
948	c ........................................................
949	IF( forward ) THEN
950
951	itau = itau + 1
952	c iday = day_ini+itau/day_step
953	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
954	c
955	c IF(time.GT.1.) THEN
956	c time = time-1.
957	c iday = iday+1
958	c ENDIF
959
960	forward = .FALSE.
961	IF( itau. EQ. itaufinp1 ) then
962	abort_message = 'Simulation finished'
963	call abort_gcm(modname,abort_message,0)
964	ENDIF
965	GO TO 2
966
967	ELSE ! of IF(forward) i.e. backward step
968
969	if (ok_iso_verif) then
970	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
971	endif !if (ok_iso_verif) then
972
973	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
974	IF(itau.EQ.itaufin) THEN
975	iav=1
976	ELSE
977	iav=0
978	ENDIF
979
980	! ! Ehouarn: re-compute geopotential for outputs
981	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
982
983	IF (ok_dynzon) THEN
984	#ifdef CPP_IOIPSL
985	c les traceurs ne sont pas sortis, trop lourd.
986	c Peut changer eventuellement si besoin.
987	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
988	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
989	& du,dudis,dutop,dufi)
990	#endif
991	ENDIF
992	IF (ok_dyn_ave) THEN
993	#ifdef CPP_IOIPSL
994	CALL writedynav(itau,vcov,
995	& ucov,teta,pk,phi,q,masse,ps,phis)
996	#endif
997	ENDIF
998
999	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1000
1001	IF(MOD(itau,iecri ).EQ.0) THEN
1002	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
1003	! ADAPTATION GCM POUR CP(T)
1004	call tpot2t(ijp1llm,teta,temp,pk)
1005	tsurpk = cpp*temp/pk
1006	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
1007	unat=0.
1008	do l=1,llm
1009	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
1010	vnat(:,l)=vcov(:,l)/cv(:)
1011	enddo
1012	#ifdef CPP_IOIPSL
1013	if (ok_dyn_ins) then
1014	! write(lunout,*) "leapfrog: call writehist (b)",
1015	! & itau,iecri
1016	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1017	endif ! of if (ok_dyn_ins)
1018	#endif
1019	! For some Grads outputs
1020	if (output_grads_dyn) then
1021	#include "write_grads_dyn.h"
1022	endif
1023
1024	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
1025
1026	c Determine whether to write to the restart.nc file
1027	c Decision can't be made in one IF statement as if
1028	c ecritstart==0 there will be a divide-by-zero error
1029	lrestart = .false.
1030	IF (itau.EQ.itaufin) THEN
1031	lrestart = .true.
1032	ELSE IF (ecritstart.GT.0) THEN
1033	IF (MOD(itau,ecritstart).EQ.0) lrestart = .true.
1034	ENDIF
1035
1036	c Write to the restart.nc file
1037	IF (lrestart) THEN
1038	if (planet_type=="mars") then
1039	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
1040	& vcov,ucov,teta,q,masse,ps)
1041	else
1042	CALL dynredem1("restart.nc",start_time,
1043	& vcov,ucov,teta,q,masse,ps)
1044	endif
1045	ENDIF ! of IF (lrestart)
1046
1047	forward = .TRUE.
1048	GO TO 1
1049
1050	ENDIF ! of IF (forward)
1051
1052	END IF ! of IF(.not.purmats)
1053
1054	STOP
1055	END

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

Download in other formats:

Original Format