Context Navigation

leapfrog_nogcm.F @ 3546

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

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

Download in other formats:

Original Format