Context Navigation

leapfrog_nogcm_pluto.F @ 3945

Last change on this file since 3945 was 3941, checked in by tbertrand, 2 months ago
PLUTO PCM : Some options in newstart for Triton simulations TB
File size: 26.6 KB

Line
1	!
2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
3	!
4	c
5	c
6	SUBROUTINE leapfrog_nogcm_pluto(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 control_mod, ONLY: planet_type,nday,day_step,iperiod,iphysiq,
17	& less1day,fractday,ndynstep,iconser,
18	& dissip_period,offline,ip_ebil_dyn,
19	& ok_dynzon,periodav,ok_dyn_ave,iecri,
20	& ok_dyn_ins,output_grads_dyn
21	use exner_hyb_m, only: exner_hyb
22	use exner_milieu_m, only: exner_milieu
23	use cpdet_mod, only: cpdet,tpot2t,t2tpot
24	use comuforc_h
25	USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs,
26	& aps,bps,presnivs,pseudoalt,preff,scaleheight
27	USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,
28	. cpp,ihf,iflag_top_bound,pi,kappa,r
29	USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
30	. statcl,conser,purmats,tidal,ok_strato
31	USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
32	. start_time,dt
33	USE callkeys_mod
34	USE tracer_h, ONLY: igcm_n2,igcm_ch4_gas,igcm_co_gas,
35	. igcm_prec_haze,igcm_haze
36
37
38	IMPLICIT NONE
39
40	c ...... Version du 30/11/24 ..........
41
42	c=======================================================================
43	c
44	c Auteur: T. Bertrand, F. Forget
45	c -------
46	c
47	c Objet:
48	c ------
49	c
50	c GCM LMD sans dynamique, pour modele saisonnier de surface
51	c
52	c=======================================================================
53	c
54	c-----------------------------------------------------------------------
55	c Declarations:
56	c -------------
57
58	#include "dimensions.h"
59	#include "paramet.h"
60	#include "comdissnew.h"
61	#include "comgeom.h"
62	!#include "com_io_dyn.h"
63	#include "iniprint.h"
64	#include "academic.h"
65
66	REAL,INTENT(IN) :: time_0 ! not used
67
68	c dynamical variables:
69	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
70	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
71	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
72	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
73	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
74	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
75	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
76
77	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
78	REAL pks(ip1jmp1) ! exner at the surface
79	REAL pk(ip1jmp1,llm) ! exner at mid-layer
80	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
81	REAL phi(ip1jmp1,llm) ! geopotential
82	REAL w(ip1jmp1,llm) ! vertical velocity
83	REAL kpd(ip1jmp1) ! TB15 = exp (-z/H)
84
85	! ADAPTATION GCM POUR CP(T)
86	REAL temp(ip1jmp1,llm) ! temperature
87	REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk
88
89	! nogcm
90	REAL tau_ps
91	real globaverage2d
92
93	c variables dynamiques intermediaire pour le transport
94	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
95
96	c variables dynamiques au pas -1
97	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
98	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
99	REAL massem1(ip1jmp1,llm)
100
101	c tendances dynamiques en */s
102	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
103	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot)
104
105	c tendances physiques */s
106	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
107	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
108
109	c variables pour le fichier histoire
110	! REAL dtav ! intervalle de temps elementaire
111
112	REAL tppn(iim),tpps(iim),tpn,tps
113	c
114	INTEGER itau,itaufinp1,iav
115	! INTEGER iday ! jour julien
116	REAL time
117
118	REAL SSUM
119
120	cym LOGICAL lafin
121	LOGICAL :: lafin=.false.
122	INTEGER ij,iq,l
123
124	REAL rdaym_ini
125	! jD_cur: jour julien courant
126	! jH_cur: heure julienne courante
127	REAL :: jD_cur, jH_cur
128	! INTEGER :: an, mois, jour
129	! REAL :: secondes
130
131	LOGICAL first,callinigrads
132	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
133	save first
134	data first/.true./
135	logical ok_sync
136	parameter (ok_sync = .true.)
137	logical physics
138
139	data callinigrads/.true./
140	character*10 string10
141
142	! REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
143	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
144
145	REAL psmean ! pression moyenne
146	REAL pqn2mean ! moyenne globale ps*qco
147	REAL pqch4mean ! moyenne globale ps*qch4
148	REAL pqcomean ! moyenne globale ps*qco
149	REAL pqhazemean
150	REAL pqprechazemean
151	REAL qmean_ch4 ! mass mean mixing ratio vap ch4
152	REAL qmean_co ! mass mean mixing ratio vap co
153	REAL qmean_haze
154	REAL qmean_prechaze
155	REAL pqch4(ip1jmp1) ! average methane mass index : ps*qch4
156	REAL pqco(ip1jmp1) ! average CO mass index : ps*q_co
157	REAL pqhaze(ip1jmp1)
158	REAL pqprechaze(ip1jmp1)
159
160	REAL p0 ! pression de reference
161	REAL p00d ! globalaverage(kpd)
162	REAL qmean_n2,qmean_n2_vert ! mass mean mixing ratio vap n2
163	REAL pqn2(ip1jmp1) ! average n2 mass index : ps*q_n2
164	REAL oldps(ip1jmp1) ! saving old pressure ps to calculate qch4
165
166	c+jld variables test conservation energie
167	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
168	C Tendance de la temp. potentiel d (theta)/ d t due a la
169	C tansformation d'energie cinetique en energie thermique
170	C cree par la dissipation
171	REAL dtetaecdt(ip1jmp1,llm)
172	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
173	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
174	CHARACTER*15 ztit
175
176	character(len=*),parameter :: modname="leapfrog"
177	character*80 abort_message
178
179	logical dissip_conservative
180	save dissip_conservative
181	data dissip_conservative/.true./
182
183	INTEGER testita
184	PARAMETER (testita = 9)
185
186	logical , parameter :: flag_verif = .false.
187
188	! for CP(T)
189	real :: dtec
190	real :: ztetaec(ip1jmp1,llm)
191
192	! TEMP : diagnostic mass
193	real :: n2mass(iip1,jjp1)
194	real :: n2ice_ij(iip1,jjp1)
195	integer :: i,j,ig
196	integer, parameter :: ngrid = 2+(jjm-1)*iim
197	! local mass
198	!real :: mq(ip1jmp1,llm)
199
200	if (nday>=0) then
201	itaufin = nday*day_step
202	else
203	! to run a given (-nday) number of dynamical steps
204	itaufin = -nday
205	endif
206	if (less1day) then
207	c MODIF VENUS: to run less than one day:
208	itaufin = int(fractday*day_step)
209	endif
210	if (ndynstep.gt.0) then
211	! running a given number of dynamical steps
212	itaufin=ndynstep
213	endif
214	itaufinp1 = itaufin +1
215
216
217
218
219	c INITIALISATIONS
220	dufi(:,:) =0.
221	dvfi(:,:) =0.
222	dtetafi(:,:)=0.
223	dqfi(:,:,:) =0.
224	dpfi(:) =0.
225
226	itau = 0
227	physics=.true.
228	if (iflag_phys==0.or.iflag_phys==2) physics=.false.
229
230	! ED18 nogcm
231	! firstcall_globaverage2d = .true.
232
233	c iday = day_ini+itau/day_step
234	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
235	c IF(time.GT.1.) THEN
236	c time = time-1.
237	c iday = iday+1
238	c ENDIF
239
240
241	c-----------------------------------------------------------------------
242	c On initialise la pression et la fonction d'Exner :
243	c --------------------------------------------------
244
245	dq(:,:,:)=0.
246	CALL pression ( ip1jmp1, ap, bp, ps, p )
247	if (pressure_exner) then
248	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
249	else
250	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
251	endif
252
253	c------------------
254	c TEST PK MONOTONE
255	c------------------
256	write(,) "Test PK"
257	do ij=1,ip1jmp1
258	do l=2,llm
259	!PRINT*,'pk(ij,l) = ',l,pk(ij,l)
260	!PRINT*,'pk(ij,l-1) = ',pk(ij,l-1)
261	if(pk(ij,l).gt.pk(ij,l-1)) then
262	c write(,) ij,l,pk(ij,l)
263	abort_message = 'PK non strictement decroissante'
264	call abort_gcm(modname,abort_message,1)
265	c write(,) "ATTENTION, Test PK deconnecte..."
266	endif
267	enddo
268	enddo
269	write(,) "Fin Test PK"
270	c stop
271
272	c------------------
273	c Preparing mixing of pressure and tracers in nogcm
274	c kpd=exp(-z/H) Needed to define a reference pressure :
275	c P0=pmean/globalaverage(kpd)
276	c thus P(i) = P0exp(-z(i)/H) = P0kpd(i)
277	c It is checked that globalaverage2d(Pi)=pmean
278	DO ij=1,ip1jmp1
279	kpd(ij) = exp(-phis(ij)/(r*38.))
280	ENDDO
281	p00d=globaverage2d(kpd) ! mean pres at ref level
282	tau_ps = tau_n2 ! constante de rappel for pressure (s)
283
284	PRINT*,'igcm_n2 = ',igcm_n2
285	do iq=1,nqtot
286	if (tname(iq)=="n2") then
287	igcm_n2=iq
288	exit
289	endif
290	enddo
291
292	c-----------------------------------------------------------------------
293	c Debut de l'integration temporelle:
294	c ----------------------------------
295
296	1 CONTINUE ! Matsuno Forward step begins here
297
298	c date: (NB: date remains unchanged for Backward step)
299	c -----
300
301	jD_cur = jD_ref + day_ini - day_ref +
302	& (itau+1)/day_step
303	jH_cur = jH_ref + start_time +
304	& mod(itau+1,day_step)/float(day_step)
305	jD_cur = jD_cur + int(jH_cur)
306	jH_cur = jH_cur - int(jH_cur)
307
308	! Save fields obtained at previous time step as '...m1'
309	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
310	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
311	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
312	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
313	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
314
315	forward = .TRUE.
316	leapf = .FALSE.
317	dt = dtvr
318
319	2 CONTINUE ! Matsuno backward or leapfrog step begins here
320
321	c gestion des appels de la physique et des dissipations:
322	c ------------------------------------------------------
323
324	apphys = .FALSE.
325	statcl = .FALSE.
326	conser = .FALSE.
327
328	IF( purmats ) THEN
329	! Purely Matsuno time stepping
330	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
331
332	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
333	s .and. physics ) apphys = .TRUE.
334	ELSE
335	! Leapfrog/Matsuno time stepping
336	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
337	IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
338	END IF
339
340	c-----------------------------------------------------------------------
341	c calcul des tendances dynamiques:
342	c --------------------------------
343
344	dv(:,:) = 0.
345	du(:,:) = 0.
346	dteta(:,:) = 0.
347	dq(:,:,:) = 0.
348
349	c-----------------------------------------------------------------------
350	c calcul des tendances physiques:
351	c -------------------------------
352	c
353	IF( purmats ) THEN
354	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
355	ELSE
356	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
357	ENDIF
358
359	IF( apphys ) THEN
360
361	CALL pression ( ip1jmp1, ap, bp, ps, p )
362	if (pressure_exner) then
363	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
364	else
365	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
366	endif
367
368	! Compute geopotential (physics might need it)
369	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
370
371	jD_cur = jD_ref + day_ini - day_ref +
372	& (itau+1)/day_step
373
374	! AS: we make jD_cur to be pday
375	jD_cur = int(day_ini + itau/day_step)
376
377	jH_cur = jH_ref + start_time +
378	& mod(itau+1,day_step)/float(day_step)
379	jH_cur = jH_ref + start_time +
380	& mod(itau,day_step)/float(day_step)
381	jD_cur = jD_cur + int(jH_cur)
382	jH_cur = jH_cur - int(jH_cur)
383
384	c Interface avec les routines de phylmd (phymars ... )
385	c -----------------------------------------------------
386
387	CALL calfis( lafin , jD_cur, jH_cur,
388	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
389	$ du,dv,dteta,dq,
390	$ flxw,
391	$ dufi,dvfi,dtetafi,dqfi,dpfi )
392
393	c ajout des tendances physiques
394	c ------------------------------
395
396	c CALL addfi( dtphys, leapf, forward ,
397	c $ ucov, vcov, teta , q ,ps ,
398	c $ dufi, dvfi, dtetafi , dqfi ,dpfi )
399
400	c CALL pression (ip1jmp1,ap,bp,ps,p)
401	c CALL massdair(p,masse)
402
403	! 1) Saving pressure and mass of N2 tracer in each mess (kg) before Horizontal mixing of pressure
404	DO ij=1,ip1jmp1
405	oldps(ij)=ps(ij)
406	ENDDO
407
408	!DO l=1, llm
409	! DO ij=1,ip1jmp1
410	! mq(ij,l) = masse(ij,l)*q(ij,l,igcm_n2)
411	! ENDDO
412	!ENDDO
413
414	! 2) Increment q with physical tendancy (mixing ratio kg/kg)
415	DO ij=1,ip1jmp1
416	if (methane) then
417	q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)+
418	& dqfi(ij,1,igcm_ch4_gas)*dtphys
419	endif
420	if (carbox) then
421	q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)+
422	& dqfi(ij,1,igcm_co_gas)*dtphys
423	endif
424	if (fasthaze) then
425	q(ij,1,igcm_haze)=q(ij,1,igcm_haze)+
426	& dqfi(ij,1,igcm_haze)*dtphys
427	q(ij,1,igcm_prec_haze)=q(ij,1,igcm_prec_haze)+
428	& dqfi(ij,1,igcm_prec_haze)*dtphys
429	endif
430	ENDDO
431
432	! 3) Update pressure
433	DO ij=1,ip1jmp1
434	ps(ij) = ps(ij) + dpfi(ij)*dtphys
435	ENDDO
436
437	c ! 4) Horizontal mixing of pressure
438	! Rappel newtonien vers psmean
439	psmean= globaverage2d(ps) ! mean pressure
440	p0=psmean/p00d
441	DO ij=1,ip1jmp1
442	ps(ij)=ps(ij) +(p0*kpd(ij)
443	& -ps(ij))*(1-exp(-dtphys/tau_ps))
444	ENDDO
445
446	write(72,*) 'psmean = ',psmean ! mean pressure redistributed
447
448	! 5) Security check: test pressure negative
449	DO ij=1,ip1jmp1
450	IF (ps(ij).le.0.) then
451	!write(,) 'Negative pressure :'
452	!write(,) 'ps = ',ps(ij),' ig = ',ij
453	!write(,) 'pmean = ',p0*kpd(ij)
454	!write(,) 'tau_ps = ',tau_ps
455	!STOP
456	ps(ij)=1.e-6*kpd(ij)/p00d
457	ENDIF
458	ENDDO
459
460	! 6) Horizontal redistribution of other tracers
461	! intermediaire de calcul si methane or CO
462	if ((methane).or.(carbox)) then
463	psmean= globaverage2d(ps)
464	endif
465
466	! 6a) Update mmr q for new pressure (both q and dqfi were calculated with old pressure in physiq)
467	DO ij=1,ip1jmp1
468	if (methane) then
469	q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)*oldps(ij)/
470	& ps(ij)
471	! average methane mass index : ps * qch4 => 'pressure of ch4'
472	pqch4(ij)= ps(ij) * q(ij,1,igcm_ch4_gas)
473	endif
474	if (carbox) then
475	q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)*oldps(ij)/
476	& ps(ij)
477	pqco(ij)= ps(ij) * q(ij,1,igcm_co_gas)
478	endif
479	if (fasthaze) then
480	q(ij,1,igcm_haze)=q(ij,1,igcm_haze)*oldps(ij)/
481	& ps(ij)
482	pqhaze(ij)= ps(ij) * q(ij,1,igcm_haze)
483	q(ij,1,igcm_prec_haze)=q(ij,1,igcm_prec_haze)*oldps(ij)/
484	& ps(ij)
485	pqprechaze(ij)= ps(ij) * q(ij,1,igcm_prec_haze)
486	endif
487
488	ENDDO
489
490	! 6b) Rappel newtonien vers q_mean
491	if (methane) then
492	pqch4mean=globaverage2d(pqch4)
493	qmean_ch4= pqch4mean / psmean
494	endif
495	if (carbox) then
496	pqcomean=globaverage2d(pqco)
497	qmean_co= pqcomean / psmean
498	endif
499	if (fasthaze) then
500	pqprechazemean=globaverage2d(pqprechaze)
501	qmean_prechaze= pqprechazemean / psmean
502	!pqhazemean=globaverage2d(pqhaze)
503	!qmean_haze= pqhazemean / psmean
504	endif
505
506	DO ij=1,ip1jmp1
507	if (methane) then
508	q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)+
509	& (qmean_ch4-q(ij,1,igcm_ch4_gas))*
510	& (1.-exp(-dtphys/tau_ch4))
511	endif
512	if (carbox) then
513	q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)+
514	& (qmean_co-q(ij,1,igcm_co_gas))*
515	& (1.-exp(-dtphys/tau_co))
516	endif ! CO
517	if (fasthaze) then
518	! q(ij,1,igcm_haze)=q(ij,1,igcm_haze)+
519	& ! (qmean_haze-q(ij,1,igcm_haze))*
520	& ! (1.-exp(-dtphys/tau_haze))
521	q(ij,1,igcm_prec_haze)=q(ij,1,igcm_prec_haze)+
522	& (qmean_prechaze-q(ij,1,igcm_prec_haze))*
523	& (1.-exp(-dtphys/tau_prechaze))
524	endif ! fasthaze
525	ENDDO
526
527	! 7) Update pressure p and pk
528	CALL pression (ip1jmp1,ap,bp,ps,p)
529	CALL massdair(p,masse)
530	if (pressure_exner) then
531	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
532	else
533	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
534	endif
535
536	ENDIF ! of IF( apphys )
537
538
539	c ********************************************************************
540	c ********************************************************************
541	c .... fin de l'integration dynamique et physique pour le pas itau ..
542	c ********************************************************************
543	c ********************************************************************
544
545	IF ( .NOT.purmats ) THEN
546	c ........................................................
547	c .............. schema matsuno + leapfrog ..............
548	c ........................................................
549
550	IF(forward. OR. leapf) THEN
551	itau= itau + 1
552	c iday= day_ini+itau/day_step
553	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
554	c IF(time.GT.1.) THEN
555	c time = time-1.
556	c iday = iday+1
557	c ENDIF
558	ENDIF
559
560	IF( itau. EQ. itaufinp1 ) then
561	if (flag_verif) then
562	write(79,*) 'ucov',ucov
563	write(80,*) 'vcov',vcov
564	write(81,*) 'teta',teta
565	write(82,*) 'ps',ps
566	write(83,*) 'q',q
567	WRITE(85,*) 'q1 = ',q(:,:,1)
568	WRITE(86,*) 'q3 = ',q(:,:,3)
569	endif
570
571	abort_message = 'Simulation finished'
572
573	call abort_gcm(modname,abort_message,0)
574	ENDIF
575
576	c-----------------------------------------------------------------------
577	c ecriture du fichier histoire moyenne:
578	c -------------------------------------
579
580	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
581	IF(itau.EQ.itaufin) THEN
582	iav=1
583	ELSE
584	iav=0
585	ENDIF
586
587	! ! Ehouarn: re-compute geopotential for outputs
588	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
589
590	IF (ok_dyn_ave) THEN
591	#ifdef CPP_IOIPSL
592	CALL writedynav(itau,vcov,
593	& ucov,teta,pk,phi,q,masse,ps,phis)
594	#endif
595	ENDIF
596
597	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
598
599	if (ok_iso_verif) then
600	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
601	endif !if (ok_iso_verif) then
602
603	c-----------------------------------------------------------------------
604	c ecriture de la bande histoire:
605	c ------------------------------
606
607	IF( MOD(itau,iecri).EQ.0) THEN
608	! Ehouarn: output only during LF or Backward Matsuno
609	if (leapf.or.(.not.leapf.and.(.not.forward))) then
610	! ADAPTATION GCM POUR CP(T)
611	call tpot2t(ijp1llm,teta,temp,pk)
612	tsurpk = cpp*temp/pk
613	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
614	unat=0.
615	do l=1,llm
616	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
617	vnat(:,l)=vcov(:,l)/cv(:)
618	enddo
619	#ifdef CPP_IOIPSL
620	if (ok_dyn_ins) then
621	! write(lunout,*) "leapfrog: call writehist, itau=",itau
622	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
623	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
624	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
625	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
626	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
627	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
628	endif ! of if (ok_dyn_ins)
629	#endif
630	! For some Grads outputs of fields
631	if (output_grads_dyn) then
632	#include "write_grads_dyn.h"
633	endif
634	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
635	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
636
637	IF(itau.EQ.itaufin) THEN
638
639	CALL dynredem1("restart.nc",start_time,
640	& vcov,ucov,teta,q,masse,ps)
641	CLOSE(99)
642	!!! Ehouarn: Why not stop here and now?
643	ENDIF ! of IF (itau.EQ.itaufin)
644
645	c-----------------------------------------------------------------------
646	c gestion de l'integration temporelle:
647	c ------------------------------------
648
649	IF( MOD(itau,iperiod).EQ.0 ) THEN
650	GO TO 1
651	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
652
653	IF( forward ) THEN
654	c fin du pas forward et debut du pas backward
655
656	forward = .FALSE.
657	leapf = .FALSE.
658	GO TO 2
659
660	ELSE
661	c fin du pas backward et debut du premier pas leapfrog
662
663	leapf = .TRUE.
664	dt = 2.*dtvr
665	GO TO 2
666	END IF ! of IF (forward)
667	ELSE
668
669	c ...... pas leapfrog .....
670
671	leapf = .TRUE.
672	dt = 2.*dtvr
673	GO TO 2
674	END IF ! of IF (MOD(itau,iperiod).EQ.0)
675	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
676
677	ELSE ! of IF (.not.purmats)
678
679	if (ok_iso_verif) then
680	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
681	endif !if (ok_iso_verif) then
682
683	c ........................................................
684	c .............. schema matsuno ...............
685	c ........................................................
686	IF( forward ) THEN
687
688	itau = itau + 1
689	c iday = day_ini+itau/day_step
690	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
691	c
692	c IF(time.GT.1.) THEN
693	c time = time-1.
694	c iday = iday+1
695	c ENDIF
696
697	forward = .FALSE.
698	IF( itau. EQ. itaufinp1 ) then
699	abort_message = 'Simulation finished'
700	call abort_gcm(modname,abort_message,0)
701	ENDIF
702	GO TO 2
703
704	ELSE ! of IF(forward) i.e. backward step
705
706	if (ok_iso_verif) then
707	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
708	endif !if (ok_iso_verif) then
709
710	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
711	IF(itau.EQ.itaufin) THEN
712	iav=1
713	ELSE
714	iav=0
715	ENDIF
716
717	! ! Ehouarn: re-compute geopotential for outputs
718	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
719
720	IF (ok_dyn_ave) THEN
721	#ifdef CPP_IOIPSL
722	CALL writedynav(itau,vcov,
723	& ucov,teta,pk,phi,q,masse,ps,phis)
724	#endif
725	ENDIF
726
727	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
728
729	IF(MOD(itau,iecri ).EQ.0) THEN
730	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
731	! ADAPTATION GCM POUR CP(T)
732	call tpot2t(ijp1llm,teta,temp,pk)
733	tsurpk = cpp*temp/pk
734	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
735	unat=0.
736	do l=1,llm
737	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
738	vnat(:,l)=vcov(:,l)/cv(:)
739	enddo
740	#ifdef CPP_IOIPSL
741	if (ok_dyn_ins) then
742	! write(lunout,*) "leapfrog: call writehist (b)",
743	! & itau,iecri
744	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
745	endif ! of if (ok_dyn_ins)
746	#endif
747	! For some Grads outputs
748	if (output_grads_dyn) then
749	#include "write_grads_dyn.h"
750	endif
751
752	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
753
754	IF(itau.EQ.itaufin) THEN
755	CALL dynredem1("restart.nc",start_time,
756	& vcov,ucov,teta,q,masse,ps)
757	ENDIF ! of IF(itau.EQ.itaufin)
758
759	forward = .TRUE.
760	GO TO 1
761
762	ENDIF ! of IF (forward)
763
764	END IF ! of IF(.not.purmats)
765
766	STOP
767	END
768
769	! ************************************************************
770	! globalaverage VERSION PLuto
771	! ************************************************************
772	FUNCTION globaverage2d(var)
773	! ************************************************************
774	IMPLICIT NONE
775	#include "dimensions.h"
776	#include "paramet.h"
777	#include "comgeom2.h"
778	!ip1jmp1 called in paramet.h = iip1 x jjp1
779	REAL var(iip1,jjp1), globaverage2d
780	integer i,j
781	REAL airetot
782	save airetot
783	logical firstcall
784	data firstcall/.true./
785	save firstcall
786
787	if (firstcall) then
788	firstcall=.false.
789	airetot =0.
790	DO j=2,jjp1-1 ! lat
791	DO i = 1, iim ! lon
792	airetot = airetot + aire(i,j)
793	END DO
794	END DO
795	DO i=1,iim
796	airetot=airetot + aire(i,1)+aire(i,jjp1)
797	ENDDO
798	end if
799
800	globaverage2d = 0.
801	DO j=2,jjp1-1
802	DO i = 1, iim
803	globaverage2d = globaverage2d + var(i,j)*aire(i,j)
804	END DO
805	END DO
806
807	DO i=1,iim
808	globaverage2d = globaverage2d + var(i,1)*aire(i,1)
809	globaverage2d = globaverage2d + var(i,jjp1)*aire(i,jjp1)
810	ENDDO
811
812	globaverage2d = globaverage2d/airetot
813	return
814	end
815

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

Download in other formats:

Original Format