Context Navigation

leapfrog.F @ 1422

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

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

Download in other formats:

Original Format