Context Navigation

leapfrog_nogcm.F @ 3521

Last change on this file since 3521 was 3422, checked in by afalco, 15 months ago
Pluto PCM: no gcm works in parallel. example restart.def for old-to-new pluto start files. AF
File size: 33.5 KB

Line
1	!
2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
3	!
4	c
5	c
6	SUBROUTINE leapfrog_nogcm(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,tname
15	USE write_field, ONLY: writefield
16	USE callkeys_mod, only: tau_n2, tau_ch4, tau_co
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 comuforc_h
26	USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs,
27	& aps,bps,presnivs,pseudoalt,preff,scaleheight
28	USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,
29	. cpp,ihf,iflag_top_bound,pi,kappa,r
30	USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
31	. statcl,conser,purmats,tidal,ok_strato
32	USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
33	. start_time,dt
34
35
36
37
38	IMPLICIT NONE
39
40	c ...... Version du 10/01/98 ..........
41
42	c avec coordonnees verticales hybrides
43	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
44
45	c=======================================================================
46	c
47	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
48	c -------
49	c
50	c Objet:
51	c ------
52	c
53	c GCM LMD nouvelle grille
54	c
55	c=======================================================================
56	c
57	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
58	c et possibilite d'appeler une fonction f(y) a derivee tangente
59	c hyperbolique a la place de la fonction a derivee sinusoidale.
60
61	c ... Possibilite de choisir le shema pour l'advection de
62	c q , en modifiant iadv dans traceur.def (10/02) .
63	c
64	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
65	c Pour Van-Leer iadv=10
66	c
67	c-----------------------------------------------------------------------
68	c Declarations:
69	c -------------
70
71	#include "dimensions.h"
72	#include "paramet.h"
73	#include "comdissnew.h"
74	#include "comgeom.h"
75	!#include "com_io_dyn.h"
76	#include "iniprint.h"
77	#include "academic.h"
78
79	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
80	! #include "clesphys.h"
81
82	REAL,INTENT(IN) :: time_0 ! not used
83
84	c dynamical variables:
85	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
86	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
87	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
88	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
89	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
90	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
91	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
92
93	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
94	REAL pks(ip1jmp1) ! exner at the surface
95	REAL pk(ip1jmp1,llm) ! exner at mid-layer
96	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
97	REAL phi(ip1jmp1,llm) ! geopotential
98	REAL w(ip1jmp1,llm) ! vertical velocity
99	REAL kpd(ip1jmp1) ! TB15 = exp (-z/H)
100
101	! ADAPTATION GCM POUR CP(T)
102	REAL temp(ip1jmp1,llm) ! temperature
103	REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk
104
105	! real zqmin,zqmax
106
107
108	! ED18 nogcm
109	REAL tau_ps
110	REAL tau_x
111	REAL tau_def
112	REAL tau_teta
113	REAL tetadpmean
114	REAL tetadp(ip1jmp1,llm)
115	REAL dpmean
116	REAL dp(ip1jmp1,llm)
117	REAL tetamean
118	real globaverage2d
119	! LOGICAL firstcall_globaverage2d
120
121
122	c variables dynamiques intermediaire pour le transport
123	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
124
125	c variables dynamiques au pas -1
126	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
127	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
128	REAL massem1(ip1jmp1,llm)
129
130	c tendances dynamiques en */s
131	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
132	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot)
133
134	c tendances de la dissipation en */s
135	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
136	REAL dtetadis(ip1jmp1,llm)
137
138	c tendances de la couche superieure */s
139	c REAL dvtop(ip1jm,llm)
140	REAL dutop(ip1jmp1,llm)
141	c REAL dtetatop(ip1jmp1,llm)
142	c REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
143
144	c TITAN : tendances due au forces de marees */s
145	REAL dvtidal(ip1jm,llm),dutidal(ip1jmp1,llm)
146
147	c tendances physiques */s
148	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
149	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
150
151	c variables pour le fichier histoire
152	! REAL dtav ! intervalle de temps elementaire
153
154	REAL tppn(iim),tpps(iim),tpn,tps
155	c
156	INTEGER itau,itaufinp1,iav
157	! INTEGER iday ! jour julien
158	REAL time
159
160	REAL SSUM
161	! REAL finvmaold(ip1jmp1,llm)
162
163	cym LOGICAL lafin
164	LOGICAL :: lafin=.false.
165	INTEGER ij,iq,l
166	! INTEGER ik
167
168	! real time_step, t_wrt, t_ops
169
170	REAL rdaym_ini
171	! jD_cur: jour julien courant
172	! jH_cur: heure julienne courante
173	REAL :: jD_cur, jH_cur
174	! INTEGER :: an, mois, jour
175	! REAL :: secondes
176
177	LOGICAL first,callinigrads
178	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
179	save first
180	data first/.true./
181	! real dt_cum
182	! character*10 infile
183	! integer zan, tau0, thoriid
184	! integer nid_ctesGCM
185	! save nid_ctesGCM
186	! real degres
187	! real rlong(iip1), rlatg(jjp1)
188	! real zx_tmp_2d(iip1,jjp1)
189	! integer ndex2d(iip1*jjp1)
190	logical ok_sync
191	parameter (ok_sync = .true.)
192	logical physics
193
194	data callinigrads/.true./
195	character*10 string10
196
197	! REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
198	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
199
200	REAL psmean ! pression moyenne
201	REAL pqxmean ! moyenne globale ps*qco
202	REAL p0 ! pression de reference
203	REAL p00d ! globalaverage(kpd)
204	REAL qmean_x,qmean_x_vert ! mass mean mixing ratio vap n2
205	REAL pqx(ip1jmp1) ! average n2 mass index : ps*q_n2
206	REAL oldps(ip1jmp1) ! saving old pressure ps to calculate qch4
207
208	c+jld variables test conservation energie
209	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
210	C Tendance de la temp. potentiel d (theta)/ d t due a la
211	C tansformation d'energie cinetique en energie thermique
212	C cree par la dissipation
213	REAL dtetaecdt(ip1jmp1,llm)
214	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
215	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
216	! REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
217	CHARACTER*15 ztit
218	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
219	!IM SAVE ip_ebil_dyn
220	!IM DATA ip_ebil_dyn/0/
221	c-jld
222
223	! integer :: itau_w ! for write_paramLMDZ_dyn.h
224
225	! character*80 dynhist_file, dynhistave_file
226	character(len=*),parameter :: modname="leapfrog"
227	character*80 abort_message
228
229	logical dissip_conservative
230	save dissip_conservative
231	data dissip_conservative/.true./
232
233	INTEGER testita
234	PARAMETER (testita = 9)
235
236	logical , parameter :: flag_verif = .false.
237
238	! for CP(T)
239	real :: dtec
240	real :: ztetaec(ip1jmp1,llm)
241
242	! TEMP : diagnostic mass
243	real :: n2mass(iip1,jjp1)
244	real :: n2ice_ij(iip1,jjp1)
245	integer,save :: igcm_n2=0 ! index of N2 tracer (if any)
246	integer,save :: igcm_ch4=0 ! index of CH4 tracer (if any)
247	integer,save :: igcm_co=0 ! index of CO tracer (if any)
248	integer :: i,j,ig
249	integer, parameter :: ngrid = 2+(jjm-1)*iim
250	! local mass
251	real :: mq(ip1jmp1,llm)
252
253	if (nday>=0) then
254	itaufin = nday*day_step
255	else
256	! to run a given (-nday) number of dynamical steps
257	itaufin = -nday
258	endif
259	if (less1day) then
260	c MODIF VENUS: to run less than one day:
261	itaufin = int(fractday*day_step)
262	endif
263	if (ndynstep.gt.0) then
264	! running a given number of dynamical steps
265	itaufin=ndynstep
266	endif
267	itaufinp1 = itaufin +1
268
269
270
271
272	c INITIALISATIONS
273	dudis(:,:) =0.
274	dvdis(:,:) =0.
275	dtetadis(:,:)=0.
276	dutop(:,:) =0.
277	c dvtop(:,:) =0.
278	c dtetatop(:,:)=0.
279	c dqtop(:,:,:) =0.
280	c dptop(:) =0.
281	dufi(:,:) =0.
282	dvfi(:,:) =0.
283	dtetafi(:,:)=0.
284	dqfi(:,:,:) =0.
285	dpfi(:) =0.
286
287	itau = 0
288	physics=.true.
289	if (iflag_phys==0.or.iflag_phys==2) physics=.false.
290
291	! ED18 nogcm
292	! firstcall_globaverage2d = .true.
293
294	c iday = day_ini+itau/day_step
295	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
296	c IF(time.GT.1.) THEN
297	c time = time-1.
298	c iday = iday+1
299	c ENDIF
300
301
302	c-----------------------------------------------------------------------
303	c On initialise la pression et la fonction d'Exner :
304	c --------------------------------------------------
305
306	dq(:,:,:)=0.
307	CALL pression ( ip1jmp1, ap, bp, ps, p )
308	if (pressure_exner) then
309	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
310	else
311	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
312	endif
313
314	c------------------
315	c TEST PK MONOTONE
316	c------------------
317	write(,) "Test PK"
318	do ij=1,ip1jmp1
319	do l=2,llm
320	! PRINT*,'pk(ij,l) = ',pk(ij,l)
321	! PRINT*,'pk(ij,l-1) = ',pk(ij,l-1)
322	if(pk(ij,l).gt.pk(ij,l-1)) then
323	c write(,) ij,l,pk(ij,l)
324	abort_message = 'PK non strictement decroissante'
325	call abort_gcm(modname,abort_message,1)
326	c write(,) "ATTENTION, Test PK deconnecte..."
327	endif
328	enddo
329	enddo
330	write(,) "Fin Test PK"
331	c stop
332	c------------------
333	c Preparing mixing of pressure and tracers in nogcm
334	c kpd=exp(-z/H) Needed to define a reference pressure :
335	c P0=pmean/globalaverage(kpd)
336	c thus P(i) = P0exp(-z(i)/H) = P0kpd(i)
337	c It is checked that globalaverage2d(Pi)=pmean
338	DO ij=1,ip1jmp1
339	kpd(ij) = exp(-phis(ij)/(r*200.))
340	ENDDO
341	p00d=globaverage2d(kpd) ! mean pres at ref level
342	tau_ps = 1. ! constante de rappel for pressure (s)
343	tau_n2 = 1 ! constante de rappel for mix ratio qn2 (s)
344	tau_def = 1.E7 ! default constante de rappel for mix ratio qX (s)
345	tau_teta = 1.E7 !constante de rappel for potentiel temperature
346
347	do iq=1,nqtot
348	if (tname(iq)=="n2") then
349	igcm_n2=iq
350	! exit
351	else if (tname(iq)=="ch4_gas") then
352	igcm_ch4=iq
353	else if (tname(iq)=="co_gas") then
354	igcm_co=iq
355	endif
356	enddo
357	PRINT*,'igcm_n2 = ',igcm_n2
358	PRINT*,'igcm_ch4 = ',igcm_ch4
359	PRINT*,'igcm_co = ',igcm_co
360
361	c-----------------------------------------------------------------------
362	c Debut de l'integration temporelle:
363	c ----------------------------------
364
365	1 CONTINUE ! Matsuno Forward step begins here
366
367	c date: (NB: date remains unchanged for Backward step)
368	c -----
369
370	jD_cur = jD_ref + day_ini - day_ref + &
371	& (itau+1)/day_step
372	jH_cur = jH_ref + start_time + &
373	& mod(itau+1,day_step)/float(day_step)
374	jD_cur = jD_cur + int(jH_cur)
375	jH_cur = jH_cur - int(jH_cur)
376
377	c
378
379	! Save fields obtained at previous time step as '...m1'
380	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
381	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
382	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
383	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
384	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
385
386	forward = .TRUE.
387	leapf = .FALSE.
388	dt = dtvr
389
390	2 CONTINUE ! Matsuno backward or leapfrog step begins here
391
392	c-----------------------------------------------------------------------
393
394	c date: (NB: only leapfrog step requires recomputing date)
395	c -----
396
397	IF (leapf) THEN
398	jD_cur = jD_ref + day_ini - day_ref +
399	& (itau+1)/day_step
400	jH_cur = jH_ref + start_time +
401	& mod(itau+1,day_step)/float(day_step)
402	jD_cur = jD_cur + int(jH_cur)
403	jH_cur = jH_cur - int(jH_cur)
404	ENDIF
405
406
407	c gestion des appels de la physique et des dissipations:
408	c ------------------------------------------------------
409	c
410	c ... P.Le Van ( 6/02/95 ) ....
411
412	! ED18: suppression des mentions de la variable apdiss dans le cas
413	! 'nogcm'
414
415	apphys = .FALSE.
416	statcl = .FALSE.
417	conser = .FALSE.
418
419
420	IF( purmats ) THEN
421	! Purely Matsuno time stepping
422	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
423
424
425	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
426	s .and. physics ) apphys = .TRUE.
427	ELSE
428	! Leapfrog/Matsuno time stepping
429	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
430
431
432	IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
433	END IF
434
435	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
436	! supress dissipation step
437
438
439
440	c-----------------------------------------------------------------------
441	c calcul des tendances dynamiques:
442	c --------------------------------
443
444	! ED18: suppression de l'onglet pour le cas nogcm
445
446
447	dv(:,:) = 0.
448	du(:,:) = 0.
449	dteta(:,:) = 0.
450	dq(:,:,:) = 0.
451
452
453
454	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
455	c
456	c-----------------------------------------------------------------------
457	c calcul des tendances physiques:
458	c -------------------------------
459	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
460	c
461	IF( purmats ) THEN
462	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
463	ELSE
464	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
465	ENDIF
466	c
467	c
468	IF( apphys ) THEN
469	c
470	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
471	c
472
473	CALL pression ( ip1jmp1, ap, bp, ps, p )
474	if (pressure_exner) then
475	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
476	else
477	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
478	endif
479
480	! Compute geopotential (physics might need it)
481	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
482
483	jD_cur = jD_ref + day_ini - day_ref + &
484	& (itau+1)/day_step
485
486	! AS: we make jD_cur to be pday
487	jD_cur = int(day_ini + itau/day_step)
488
489	! print*,'itau =',itau
490	! print*,'day_step =',day_step
491	! print*,'itau/day_step =',itau/day_step
492
493
494	jH_cur = jH_ref + start_time + &
495	& mod(itau+1,day_step)/float(day_step)
496	jH_cur = jH_ref + start_time + &
497	& mod(itau,day_step)/float(day_step)
498	jD_cur = jD_cur + int(jH_cur)
499	jH_cur = jH_cur - int(jH_cur)
500
501
502
503	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
504	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
505	! write(lunout,*)'current date = ',an, mois, jour, secondes
506
507	c rajout debug
508	c lafin = .true.
509
510
511	c Interface avec les routines de phylmd (phymars ... )
512	c -----------------------------------------------------
513
514	c+jld
515
516	c Diagnostique de conservation de l'Energie : initialisation
517	IF (ip_ebil_dyn.ge.1 ) THEN
518	ztit='bil dyn'
519	! Ehouarn: be careful, diagedyn is Earth-specific!
520	IF (planet_type.eq."earth") THEN
521	CALL diagedyn(ztit,2,1,1,dtphys
522	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
523	ENDIF
524	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
525	c-jld
526	#ifdef CPP_IOIPSL
527	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
528	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
529	c IF (first) THEN
530	c first=.false.
531	c#include "ini_paramLMDZ_dyn.h"
532	c ENDIF
533	c
534	c#include "write_paramLMDZ_dyn.h"
535	c
536	#endif
537	! #endif of #ifdef CPP_IOIPSL
538
539	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
540	! print*,'---------LEAPFROG---------------'
541	! print*,''
542	! print*,'> AVANT CALFIS'
543	! print*,''
544	! print*,'teta(3049,:) = ',teta(3049,:)
545	! print*,''
546	! print*,'dteta(3049,:) = ',dteta(3049,:)
547	! print*,''
548	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
549	! print*,''
550
551	CALL calfis( lafin , jD_cur, jH_cur,
552	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
553	$ du,dv,dteta,dq,
554	$ flxw,
555	$ dufi,dvfi,dtetafi,dqfi,dpfi )
556
557	c call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
558	c call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
559	c call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
560	! print*,'> APRES CALFIS (AVANT ADDFI)'
561	! print*,''
562	! print*,'teta(3049,:) = ',teta(3049,:)
563	! print*,''
564	! print*,'dteta(3049,:) = ',dteta(3049,:)
565	! print*,''
566	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
567	! print*,''
568
569
570	c ajout des tendances physiques
571	c ------------------------------
572
573	CALL addfi( dtphys, leapf, forward ,
574	$ ucov, vcov, teta , q ,ps ,
575	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
576
577	! print*,'> APRES ADDFI'
578	! print*,''
579	! print*,'teta(3049,:) = ',teta(3049,:)
580	! print*,''
581	! print*,'dteta(3049,:) = ',dteta(3049,:)
582	! print*,''
583	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
584	! print*,''
585
586	CALL pression (ip1jmp1,ap,bp,ps,p)
587	CALL massdair(p,masse)
588
589	! Mass of tracers in each mess (kg) before Horizontal mixing of pressure
590	c -------------------------------------------------------------------
591
592	DO l=1, llm
593	DO ij=1,ip1jmp1
594	mq(ij,l) = masse(ij,l)*q(ij,l,igcm_n2)
595	ENDDO
596	ENDDO
597
598	c Horizontal mixing of pressure
599	c ------------------------------
600	! Rappel newtonien vers psmean
601	psmean= globaverage2d(ps) ! mean pressure
602	p0=psmean/p00d
603	DO ij=1,ip1jmp1
604	oldps(ij)=ps(ij)
605	ps(ij)=ps(ij) +(p0*kpd(ij)
606	& -ps(ij))*(1-exp(-dtphys/tau_ps))
607	ENDDO
608
609	write(72,*) 'psmean = ',psmean ! mean pressure redistributed
610
611	! security check: test pressure negative
612	DO ij=1,ip1jmp1
613	IF (ps(ij).le.0.) then
614	!write(,) 'Negative pressure :'
615	!write(,) 'ps = ',ps(ij),' ig = ',ij
616	!write(,) 'pmean = ',p0*kpd(ij)
617	!write(,) 'tau_ps = ',tau_ps
618	!STOP
619	ps(ij)=0.0000001*kpd(ij)/p00d
620	ENDIF
621	ENDDO
622	!***********************
623	! Correction on teta due to surface pressure changes
624	DO l = 1,llm
625	DO ij = 1,ip1jmp1
626	teta(ij,l)= teta(ij,l)(ps(ij)/oldps(ij))*kappa
627	ENDDO
628	ENDDO
629	!***********************
630
631
632
633	! ! update pressure and update p and pk
634	! DO ij=1,ip1jmp1
635	! ps(ij) = ps(ij) + dpfi(ij)*dtphys
636	! ENDDO
637	CALL pression (ip1jmp1,ap,bp,ps,p)
638	CALL massdair(p,masse)
639	if (pressure_exner) then
640	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
641	else
642	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
643	endif
644
645	! Update tracers after Horizontal mixing of pressure to ! conserve tracer mass
646	c -------------------------------------------------------------------
647	DO l=1, llm
648	DO ij=1,ip1jmp1
649	q(ij,l,igcm_n2) = mq(ij,l)/ masse(ij,l)
650	ENDDO
651	ENDDO
652
653
654
655	c Horizontal mixing of pressure
656	c ------------------------------
657	! Rappel newtonien vers psmean
658	psmean= globaverage2d(ps) ! mean pressure
659	! ! increment q_n2 with physical tendancy
660	! IF (igcm_n2.ne.0) then
661	! DO l=1, llm
662	! DO ij=1,ip1jmp1
663	! q(ij,l,igcm_n2)=q(ij,l,igcm_n2)+
664	! & dqfi(ij,l,igcm_n2)*dtphys
665	! ENDDO
666	! ENDDO
667	! ENDIF
668
669	c Mixing N2 vertically ! not used for pluto ?
670	c --------------------------
671	! if (igcm_n2.ne.0) then
672	! DO ij=1,ip1jmp1
673	! qmean_x_vert=0.
674	! DO l=1, llm
675	! qmean_x_vert= qmean_x_vert
676	! & + q(ij,l,igcm_n2)*( p(ij,l) - p(ij,l+1))
677	! END DO
678	! qmean_x_vert= qmean_x_vert/ps(ij)
679	! DO l=1, llm
680	! q(ij,l,igcm_n2)= qmean_x_vert
681	! END DO
682	! END DO
683	! end if
684
685
686
687	c Horizontal mixing of pressure
688	c ------------------------------
689	! Rappel newtonien vers psmean
690	psmean= globaverage2d(ps) ! mean pressure
691
692	c Horizontal mixing tracers, and temperature (replace dynamics in nogcm)
693	c ---------------------------------------------------------------------------------
694
695	! Simulate redistribution by dynamics for qX
696	DO iq=1,nqtot
697	if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ch4).or.
698	& (iq.eq.igcm_co)) then
699
700	DO ij=1,ip1jmp1
701	pqx(ij)= ps(ij) * q(ij,1,iq)
702	ENDDO
703	pqxmean=globaverage2d(pqx)
704
705	! Rappel newtonien vers qx_mean
706	qmean_x= pqxmean / psmean
707
708	tau_x = tau_def
709	if (iq.eq.igcm_n2) then
710	tau_x = tau_n2
711	else if (iq.eq.igcm_ch4) then
712	tau_x = tau_ch4
713	else if (iq.eq.igcm_co) then
714	tau_x = tau_co
715	end if
716	PRINT*,' tau_x ',iq,tau_x
717
718	DO ij=1,ip1jmp1
719	q(ij,1,iq)=q(ij,1,iq)+
720	& (qmean_x-q(ij,1,iq))*
721	& (1.-exp(-dtphys/tau_x))
722	ENDDO
723
724	DO l=2, llm
725	DO ij=1,ip1jmp1
726	q(ij,l,iq)=q(ij,1,iq)
727	END DO
728	END DO
729
730	! TEMPORAIRE (ED)
731	! PRINT*,'psmean = ',psmean
732	! PRINT*,'qmean_x = ',qmean_x
733	! PRINT*,'pqxmean = ',pqxmean
734	! PRINT*,' q(50,1) = ',iq,q(50,1,iq)
735	! PRINT*,' q(50,2) = ',iq,q(50,2,iq)
736	! PRINT*,' q(50,3) = ',iq,q(50,3,iq)
737
738	endif ! igcm_n2.ne.0
739	enddo
740
741
742	! ***************************************s*************
743
744	c Horizontal mixing of pressure
745	c ------------------------------
746	! Rappel newtonien vers psmean
747	psmean= globaverage2d(ps) ! mean pressure
748
749
750	c Mixing Temperature horizontally
751	c -------------------------------
752	! initialize variables that will be averaged
753	DO l=1,llm
754	DO ij=1,ip1jmp1
755	dp(ij,l) = p(ij,l) - p(ij,l+1)
756	tetadp(ij,l) = teta(ij,l)*dp(ij,l)
757	ENDDO
758	ENDDO
759
760	DO l=1,llm
761	tetadpmean = globaverage2d(tetadp(:,l))
762	dpmean = globaverage2d(dp(:,l))
763	tetamean = tetadpmean / dpmean
764	DO ij=1,ip1jmp1
765	teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) *
766	& (1 - exp(-dtphys/tau_teta))
767	ENDDO
768	ENDDO
769
770
771	ENDIF ! of IF( apphys )
772
773
774	c ********************************************************************
775	c ********************************************************************
776	c .... fin de l'integration dynamique et physique pour le pas itau ..
777	c ********************************************************************
778	c ********************************************************************
779
780
781	IF ( .NOT.purmats ) THEN
782	c ........................................................
783	c .............. schema matsuno + leapfrog ..............
784	c ........................................................
785
786	IF(forward. OR. leapf) THEN
787	itau= itau + 1
788	c iday= day_ini+itau/day_step
789	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
790	c IF(time.GT.1.) THEN
791	c time = time-1.
792	c iday = iday+1
793	c ENDIF
794	ENDIF
795
796
797	IF( itau. EQ. itaufinp1 ) then
798	if (flag_verif) then
799	write(79,*) 'ucov',ucov
800	write(80,*) 'vcov',vcov
801	write(81,*) 'teta',teta
802	write(82,*) 'ps',ps
803	write(83,*) 'q',q
804	WRITE(85,*) 'q1 = ',q(:,:,1)
805	WRITE(86,*) 'q3 = ',q(:,:,3)
806	endif
807
808	abort_message = 'Simulation finished'
809
810	call abort_gcm(modname,abort_message,0)
811	ENDIF
812	c-----------------------------------------------------------------------
813	c ecriture du fichier histoire moyenne:
814	c -------------------------------------
815
816	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
817	IF(itau.EQ.itaufin) THEN
818	iav=1
819	ELSE
820	iav=0
821	ENDIF
822
823	! ! Ehouarn: re-compute geopotential for outputs
824	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
825
826	IF (ok_dynzon) THEN
827	#ifdef CPP_IOIPSL
828	c les traceurs ne sont pas sortis, trop lourd.
829	c Peut changer eventuellement si besoin.
830	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
831	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
832	& du,dudis,dutop,dufi)
833	#endif
834	END IF
835	IF (ok_dyn_ave) THEN
836	#ifdef CPP_IOIPSL
837	CALL writedynav(itau,vcov,
838	& ucov,teta,pk,phi,q,masse,ps,phis)
839	#endif
840	ENDIF
841
842	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
843
844	if (ok_iso_verif) then
845	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
846	endif !if (ok_iso_verif) then
847
848	c-----------------------------------------------------------------------
849	c ecriture de la bande histoire:
850	c ------------------------------
851
852	IF( MOD(itau,iecri).EQ.0) THEN
853	! Ehouarn: output only during LF or Backward Matsuno
854	if (leapf.or.(.not.leapf.and.(.not.forward))) then
855	! ADAPTATION GCM POUR CP(T)
856	call tpot2t(ijp1llm,teta,temp,pk)
857	tsurpk = cpp*temp/pk
858	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
859	unat=0.
860	do l=1,llm
861	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
862	vnat(:,l)=vcov(:,l)/cv(:)
863	enddo
864	#ifdef CPP_IOIPSL
865	if (ok_dyn_ins) then
866	! write(lunout,*) "leapfrog: call writehist, itau=",itau
867	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
868	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
869	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
870	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
871	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
872	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
873	endif ! of if (ok_dyn_ins)
874	#endif
875	! For some Grads outputs of fields
876	if (output_grads_dyn) then
877	#include "write_grads_dyn.h"
878	endif
879	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
880	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
881
882	IF(itau.EQ.itaufin) THEN
883
884	CALL dynredem1("restart.nc",start_time,
885	& vcov,ucov,teta,q,masse,ps)
886	CLOSE(99)
887	!!! Ehouarn: Why not stop here and now?
888	ENDIF ! of IF (itau.EQ.itaufin)
889
890	c-----------------------------------------------------------------------
891	c gestion de l'integration temporelle:
892	c ------------------------------------
893
894	IF( MOD(itau,iperiod).EQ.0 ) THEN
895	GO TO 1
896	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
897
898	IF( forward ) THEN
899	c fin du pas forward et debut du pas backward
900
901	forward = .FALSE.
902	leapf = .FALSE.
903	GO TO 2
904
905	ELSE
906	c fin du pas backward et debut du premier pas leapfrog
907
908	leapf = .TRUE.
909	dt = 2.*dtvr
910	GO TO 2
911	END IF ! of IF (forward)
912	ELSE
913
914	c ...... pas leapfrog .....
915
916	leapf = .TRUE.
917	dt = 2.*dtvr
918	GO TO 2
919	END IF ! of IF (MOD(itau,iperiod).EQ.0)
920	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
921
922	ELSE ! of IF (.not.purmats)
923
924	if (ok_iso_verif) then
925	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
926	endif !if (ok_iso_verif) then
927
928	c ........................................................
929	c .............. schema matsuno ...............
930	c ........................................................
931	IF( forward ) THEN
932
933	itau = itau + 1
934	c iday = day_ini+itau/day_step
935	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
936	c
937	c IF(time.GT.1.) THEN
938	c time = time-1.
939	c iday = iday+1
940	c ENDIF
941
942	forward = .FALSE.
943	IF( itau. EQ. itaufinp1 ) then
944	abort_message = 'Simulation finished'
945	call abort_gcm(modname,abort_message,0)
946	ENDIF
947	GO TO 2
948
949	ELSE ! of IF(forward) i.e. backward step
950
951	if (ok_iso_verif) then
952	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
953	endif !if (ok_iso_verif) then
954
955	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
956	IF(itau.EQ.itaufin) THEN
957	iav=1
958	ELSE
959	iav=0
960	ENDIF
961
962	! ! Ehouarn: re-compute geopotential for outputs
963	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
964
965	IF (ok_dynzon) THEN
966	#ifdef CPP_IOIPSL
967	c les traceurs ne sont pas sortis, trop lourd.
968	c Peut changer eventuellement si besoin.
969	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
970	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
971	& du,dudis,dutop,dufi)
972	#endif
973	ENDIF
974	IF (ok_dyn_ave) THEN
975	#ifdef CPP_IOIPSL
976	CALL writedynav(itau,vcov,
977	& ucov,teta,pk,phi,q,masse,ps,phis)
978	#endif
979	ENDIF
980
981	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
982
983	IF(MOD(itau,iecri ).EQ.0) THEN
984	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
985	! ADAPTATION GCM POUR CP(T)
986	call tpot2t(ijp1llm,teta,temp,pk)
987	tsurpk = cpp*temp/pk
988	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
989	unat=0.
990	do l=1,llm
991	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
992	vnat(:,l)=vcov(:,l)/cv(:)
993	enddo
994	#ifdef CPP_IOIPSL
995	if (ok_dyn_ins) then
996	! write(lunout,*) "leapfrog: call writehist (b)",
997	! & itau,iecri
998	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
999	endif ! of if (ok_dyn_ins)
1000	#endif
1001	! For some Grads outputs
1002	if (output_grads_dyn) then
1003	#include "write_grads_dyn.h"
1004	endif
1005
1006	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
1007
1008	IF(itau.EQ.itaufin) THEN
1009	CALL dynredem1("restart.nc",start_time,
1010	& vcov,ucov,teta,q,masse,ps)
1011	ENDIF ! of IF(itau.EQ.itaufin)
1012
1013	forward = .TRUE.
1014	GO TO 1
1015
1016	ENDIF ! of IF (forward)
1017
1018	END IF ! of IF(.not.purmats)
1019
1020	STOP
1021	END
1022
1023	! ************************************************************
1024	! globalaverage VERSION PLuto
1025	! ************************************************************
1026	FUNCTION globaverage2d(var)
1027	! ************************************************************
1028	IMPLICIT NONE
1029	#include "dimensions.h"
1030	#include "paramet.h"
1031	#include "comgeom2.h"
1032	!ip1jmp1 called in paramet.h = iip1 x jjp1
1033	REAL var(iip1,jjp1), globaverage2d
1034	integer i,j
1035	REAL airetot
1036	save airetot
1037	logical firstcall
1038	data firstcall/.true./
1039	save firstcall
1040
1041	if (firstcall) then
1042	firstcall=.false.
1043	airetot =0.
1044	DO j=2,jjp1-1 ! lat
1045	DO i = 1, iim ! lon
1046	airetot = airetot + aire(i,j)
1047	END DO
1048	END DO
1049	DO i=1,iim
1050	airetot=airetot + aire(i,1)+aire(i,jjp1)
1051	ENDDO
1052	end if
1053
1054	globaverage2d = 0.
1055	DO j=2,jjp1-1
1056	DO i = 1, iim
1057	globaverage2d = globaverage2d + var(i,j)*aire(i,j)
1058	END DO
1059	END DO
1060
1061	DO i=1,iim
1062	globaverage2d = globaverage2d + var(i,1)*aire(i,1)
1063	globaverage2d = globaverage2d + var(i,jjp1)*aire(i,jjp1)
1064	ENDDO
1065
1066	globaverage2d = globaverage2d/airetot
1067	return
1068	end
1069

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

Download in other formats:

Original Format