Context Navigation

leapfrog.F @ 1056

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

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

Download in other formats:

Original Format