Context Navigation

leapfrog.F @ 1107

Last change on this file since 1107 was 1107, checked in by emillour, 11 years ago

Common dynamics: Updates and modifications to enable running Mars physics with

LMDZ.COMMON dynamics:

For compilation: adapted makelmdz, create_make_gcm and makelmdz_fcm, bld.cfg to compile aeronomy routines in "aerono$physique" if it exists, and added "-P -traditional" preprocessing flags in "arch-linux-ifort*"
Added function "cbrt.F" (cubic root) in 'bibio'
Adapted the reading/writing of dynamics (re)start.nc files for Mars. The main issue is that different information (on time, reference and current) is stored and used differently, hence a few if (planet_type =="mars") here and there. Moreover in the martian case there is the possibility to store fields over multiple times. Some Mars-specific variables (ecritphy,ecritstart,timestart) added in control_mod.F and (hour_ini) in temps.h

File size: 29.7 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").or.
487	& (planet_type .eq."mars")) THEN
488	! AS: we make jD_cur to be pday
489	jD_cur = int(day_ini + itau/day_step)
490	ENDIF
491
492	jH_cur = jH_ref + start_time + &
493	& mod(itau,day_step)/float(day_step)
494	jD_cur = jD_cur + int(jH_cur)
495	jH_cur = jH_cur - int(jH_cur)
496	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
497	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
498	! write(lunout,*)'current date = ',an, mois, jour, secondes
499
500	c rajout debug
501	c lafin = .true.
502
503
504	c Interface avec les routines de phylmd (phymars ... )
505	c -----------------------------------------------------
506
507	c+jld
508
509	c Diagnostique de conservation de l'énergie : initialisation
510	IF (ip_ebil_dyn.ge.1 ) THEN
511	ztit='bil dyn'
512	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
513	IF (planet_type.eq."earth") THEN
514	CALL diagedyn(ztit,2,1,1,dtphys
515	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
516	ENDIF
517	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
518	c-jld
519	#ifdef CPP_IOIPSL
520	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
521	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
522	c IF (first) THEN
523	c first=.false.
524	c#include "ini_paramLMDZ_dyn.h"
525	c ENDIF
526	c
527	c#include "write_paramLMDZ_dyn.h"
528	c
529	#endif
530	! #endif of #ifdef CPP_IOIPSL
531
532	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
533
534	CALL calfis( lafin , jD_cur, jH_cur,
535	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
536	$ du,dv,dteta,dq,
537	$ flxw,
538	$ dufi,dvfi,dtetafi,dqfi,dpfi )
539
540	c call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
541	c call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
542	c call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
543
544	c ajout des tendances physiques:
545	c ------------------------------
546	CALL addfi( dtphys, leapf, forward ,
547	$ ucov, vcov, teta , q ,ps ,
548	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
549
550	c Couche superieure :
551	c -------------------
552	IF (iflag_top_bound > 0) THEN
553	CALL top_bound(vcov,ucov,teta,masse,dtphys,dutop)
554	dutop(:,:)=dutop(:,:)/dtphys ! convert to a tendency in (m/s)/s
555	ENDIF
556
557	c Diagnostique de conservation de l'énergie : difference
558	IF (ip_ebil_dyn.ge.1 ) THEN
559	ztit='bil phys'
560	IF (planet_type.eq."earth") THEN
561	CALL diagedyn(ztit,2,1,1,dtphys
562	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
563	ENDIF
564	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
565
566	ENDIF ! of IF( apphys )
567
568	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
569	! Academic case : Simple friction and Newtonan relaxation
570	! -------------------------------------------------------
571	DO l=1,llm
572	DO ij=1,ip1jmp1
573	teta(ij,l)=teta(ij,l)-dtvr*
574	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
575	ENDDO
576	ENDDO ! of DO l=1,llm
577
578	if (planet_type.eq."giant") then
579	! add an intrinsic heat flux at the base of the atmosphere
580	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
581	endif
582
583	call friction(ucov,vcov,dtvr)
584
585
586	! Sponge layer (if any)
587	IF (ok_strato) THEN
588	CALL top_bound(vcov,ucov,teta,masse,dtvr,dutop)
589	dutop(:,:)=dutop(:,:)/dtvr ! convert to a tendency in (m/s)/s
590	ENDIF ! of IF (ok_strato)
591	ENDIF ! of IF (iflag_phys.EQ.2)
592
593
594	c-jld
595
596	CALL pression ( ip1jmp1, ap, bp, ps, p )
597	if (pressure_exner) then
598	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
599	else
600	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
601	endif
602
603
604	c-----------------------------------------------------------------------
605	c dissipation horizontale et verticale des petites echelles:
606	c ----------------------------------------------------------
607
608	IF(apdiss) THEN
609
610	! sponge layer
611	if (callsponge) then
612	CALL sponge(ucov,vcov,teta,ps,dtdiss,mode_sponge)
613	endif
614
615	c calcul de l'energie cinetique avant dissipation
616	call covcont(llm,ucov,vcov,ucont,vcont)
617	call enercin(vcov,ucov,vcont,ucont,ecin0)
618
619	c dissipation
620	! ADAPTATION GCM POUR CP(T)
621	call tpot2t(ijp1llm,teta,temp,pk)
622
623	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
624	ucov=ucov+dudis
625	vcov=vcov+dvdis
626	dudis=dudis/dtdiss ! passage en (m/s)/s
627	dvdis=dvdis/dtdiss ! passage en (m/s)/s
628
629	c------------------------------------------------------------------------
630	if (dissip_conservative) then
631	C On rajoute la tendance due a la transform. Ec -> E therm. cree
632	C lors de la dissipation
633	call covcont(llm,ucov,vcov,ucont,vcont)
634	call enercin(vcov,ucov,vcont,ucont,ecin)
635	! ADAPTATION GCM POUR CP(T)
636	do ij=1,ip1jmp1
637	do l=1,llm
638	dtec = (ecin0(ij,l)-ecin(ij,l))/cpdet(temp(ij,l))
639	temp(ij,l) = temp(ij,l) + dtec
640	enddo
641	enddo
642	call t2tpot(ijp1llm,temp,ztetaec,pk)
643	dtetaecdt=ztetaec-teta
644	dtetadis=dtetadis+dtetaecdt
645	endif
646	teta=teta+dtetadis
647	dtetadis=dtetadis/dtdiss ! passage en K/s
648	c------------------------------------------------------------------------
649
650
651	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
652	c ... Calcul de la valeur moyenne, unique de h aux poles .....
653	c
654
655	DO l = 1, llm
656	DO ij = 1,iim
657	tppn(ij) = aire( ij ) * teta( ij ,l)
658	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
659	ENDDO
660	tpn = SSUM(iim,tppn,1)/apoln
661	tps = SSUM(iim,tpps,1)/apols
662
663	DO ij = 1, iip1
664	teta( ij ,l) = tpn
665	teta(ij+ip1jm,l) = tps
666	ENDDO
667	ENDDO
668
669	if (1 == 0) then
670	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
671	!!! 2) should probably not be here anyway
672	!!! but are kept for those who would want to revert to previous behaviour
673	DO ij = 1,iim
674	tppn(ij) = aire( ij ) * ps ( ij )
675	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
676	ENDDO
677	tpn = SSUM(iim,tppn,1)/apoln
678	tps = SSUM(iim,tpps,1)/apols
679
680	DO ij = 1, iip1
681	ps( ij ) = tpn
682	ps(ij+ip1jm) = tps
683	ENDDO
684	endif ! of if (1 == 0)
685
686	END IF ! of IF(apdiss)
687
688	c ajout debug
689	c IF( lafin ) then
690	c abort_message = 'Simulation finished'
691	c call abort_gcm(modname,abort_message,0)
692	c ENDIF
693
694	c ********************************************************************
695	c ********************************************************************
696	c .... fin de l'integration dynamique et physique pour le pas itau ..
697	c ********************************************************************
698	c ********************************************************************
699
700	c preparation du pas d'integration suivant ......
701
702	IF ( .NOT.purmats ) THEN
703	c ........................................................
704	c .............. schema matsuno + leapfrog ..............
705	c ........................................................
706
707	IF(forward. OR. leapf) THEN
708	itau= itau + 1
709	c iday= day_ini+itau/day_step
710	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
711	c IF(time.GT.1.) THEN
712	c time = time-1.
713	c iday = iday+1
714	c ENDIF
715	ENDIF
716
717
718	IF( itau. EQ. itaufinp1 ) then
719	if (flag_verif) then
720	write(79,*) 'ucov',ucov
721	write(80,*) 'vcov',vcov
722	write(81,*) 'teta',teta
723	write(82,*) 'ps',ps
724	write(83,*) 'q',q
725	WRITE(85,*) 'q1 = ',q(:,:,1)
726	WRITE(86,*) 'q3 = ',q(:,:,3)
727	endif
728
729	abort_message = 'Simulation finished'
730
731	call abort_gcm(modname,abort_message,0)
732	ENDIF
733	c-----------------------------------------------------------------------
734	c ecriture du fichier histoire moyenne:
735	c -------------------------------------
736
737	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
738	IF(itau.EQ.itaufin) THEN
739	iav=1
740	ELSE
741	iav=0
742	ENDIF
743
744	IF (ok_dynzon) THEN
745	#ifdef CPP_IOIPSL
746	c les traceurs ne sont pas sortis, trop lourd.
747	c Peut changer eventuellement si besoin.
748	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
749	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
750	& du,dudis,dutop,dufi)
751	#endif
752	END IF
753	IF (ok_dyn_ave) THEN
754	#ifdef CPP_IOIPSL
755	CALL writedynav(itau,vcov,
756	& ucov,teta,pk,phi,q,masse,ps,phis)
757	#endif
758	ENDIF
759
760	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
761
762	c-----------------------------------------------------------------------
763	c ecriture de la bande histoire:
764	c ------------------------------
765
766	IF( MOD(itau,iecri).EQ.0) THEN
767	! Ehouarn: output only during LF or Backward Matsuno
768	if (leapf.or.(.not.leapf.and.(.not.forward))) then
769	! ADAPTATION GCM POUR CP(T)
770	call tpot2t(ijp1llm,teta,temp,pk)
771	tsurpk = cpp*temp/pk
772	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
773	unat=0.
774	do l=1,llm
775	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
776	vnat(:,l)=vcov(:,l)/cv(:)
777	enddo
778	#ifdef CPP_IOIPSL
779	if (ok_dyn_ins) then
780	! write(lunout,*) "leapfrog: call writehist, itau=",itau
781	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
782	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
783	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
784	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
785	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
786	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
787	endif ! of if (ok_dyn_ins)
788	#endif
789	! For some Grads outputs of fields
790	if (output_grads_dyn) then
791	#include "write_grads_dyn.h"
792	endif
793	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
794	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
795
796	IF(itau.EQ.itaufin) THEN
797
798	if (planet_type=="mars") then
799	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
800	& vcov,ucov,teta,q,masse,ps)
801	else
802	CALL dynredem1("restart.nc",start_time,
803	& vcov,ucov,teta,q,masse,ps)
804	endif
805	CLOSE(99)
806	!!! Ehouarn: Why not stop here and now?
807	ENDIF ! of IF (itau.EQ.itaufin)
808
809	c-----------------------------------------------------------------------
810	c gestion de l'integration temporelle:
811	c ------------------------------------
812
813	IF( MOD(itau,iperiod).EQ.0 ) THEN
814	GO TO 1
815	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
816
817	IF( forward ) THEN
818	c fin du pas forward et debut du pas backward
819
820	forward = .FALSE.
821	leapf = .FALSE.
822	GO TO 2
823
824	ELSE
825	c fin du pas backward et debut du premier pas leapfrog
826
827	leapf = .TRUE.
828	dt = 2.*dtvr
829	GO TO 2
830	END IF ! of IF (forward)
831	ELSE
832
833	c ...... pas leapfrog .....
834
835	leapf = .TRUE.
836	dt = 2.*dtvr
837	GO TO 2
838	END IF ! of IF (MOD(itau,iperiod).EQ.0)
839	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
840
841	ELSE ! of IF (.not.purmats)
842
843	c ........................................................
844	c .............. schema matsuno ...............
845	c ........................................................
846	IF( forward ) THEN
847
848	itau = itau + 1
849	c iday = day_ini+itau/day_step
850	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
851	c
852	c IF(time.GT.1.) THEN
853	c time = time-1.
854	c iday = iday+1
855	c ENDIF
856
857	forward = .FALSE.
858	IF( itau. EQ. itaufinp1 ) then
859	abort_message = 'Simulation finished'
860	call abort_gcm(modname,abort_message,0)
861	ENDIF
862	GO TO 2
863
864	ELSE ! of IF(forward) i.e. backward step
865
866	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
867	IF(itau.EQ.itaufin) THEN
868	iav=1
869	ELSE
870	iav=0
871	ENDIF
872
873	IF (ok_dynzon) THEN
874	#ifdef CPP_IOIPSL
875	c les traceurs ne sont pas sortis, trop lourd.
876	c Peut changer eventuellement si besoin.
877	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
878	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
879	& du,dudis,dutop,dufi)
880	#endif
881	ENDIF
882	IF (ok_dyn_ave) THEN
883	#ifdef CPP_IOIPSL
884	CALL writedynav(itau,vcov,
885	& ucov,teta,pk,phi,q,masse,ps,phis)
886	#endif
887	ENDIF
888
889	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
890
891	IF(MOD(itau,iecri ).EQ.0) THEN
892	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
893	! ADAPTATION GCM POUR CP(T)
894	call tpot2t(ijp1llm,teta,temp,pk)
895	tsurpk = cpp*temp/pk
896	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
897	unat=0.
898	do l=1,llm
899	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
900	vnat(:,l)=vcov(:,l)/cv(:)
901	enddo
902	#ifdef CPP_IOIPSL
903	if (ok_dyn_ins) then
904	! write(lunout,*) "leapfrog: call writehist (b)",
905	! & itau,iecri
906	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
907	endif ! of if (ok_dyn_ins)
908	#endif
909	! For some Grads outputs
910	if (output_grads_dyn) then
911	#include "write_grads_dyn.h"
912	endif
913
914	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
915
916	IF(itau.EQ.itaufin) THEN
917	if (planet_type=="mars") then
918	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
919	& vcov,ucov,teta,q,masse,ps)
920	else
921	CALL dynredem1("restart.nc",start_time,
922	& vcov,ucov,teta,q,masse,ps)
923	endif
924	ENDIF ! of IF(itau.EQ.itaufin)
925
926	forward = .TRUE.
927	GO TO 1
928
929	ENDIF ! of IF (forward)
930
931	END IF ! of IF(.not.purmats)
932
933	STOP
934	END

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

Download in other formats:

Original Format