Context Navigation

leapfrog_nogcm.F @ 3000

Last change on this file since 3000 was 2108, checked in by emillour, 6 years ago
Follow-up of r2016 commit of nogcm main program: adapt leapfrog_gcm so that it does not explicitely depend on "phymars" and thus causes compilation issues when using other physics packages. EM
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 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	& aps,bps,presnivs,pseudoalt,preff,scaleheight
29	USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,
30	. cpp,ihf,iflag_top_bound,pi,kappa,r
31	USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
32	. statcl,conser,purmats,tidal,ok_strato
33	USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
34	. start_time,dt
35
36
37
38
39	IMPLICIT NONE
40
41	c ...... Version du 10/01/98 ..........
42
43	c avec coordonnees verticales hybrides
44	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
45
46	c=======================================================================
47	c
48	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
49	c -------
50	c
51	c Objet:
52	c ------
53	c
54	c GCM LMD nouvelle grille
55	c
56	c=======================================================================
57	c
58	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
59	c et possibilite d'appeler une fonction f(y) a derivee tangente
60	c hyperbolique a la place de la fonction a derivee sinusoidale.
61
62	c ... Possibilite de choisir le shema pour l'advection de
63	c q , en modifiant iadv dans traceur.def (10/02) .
64	c
65	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
66	c Pour Van-Leer iadv=10
67	c
68	c-----------------------------------------------------------------------
69	c Declarations:
70	c -------------
71
72	#include "dimensions.h"
73	#include "paramet.h"
74	#include "comdissnew.h"
75	#include "comgeom.h"
76	!#include "com_io_dyn.h"
77	#include "iniprint.h"
78	#include "academic.h"
79
80	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
81	! #include "clesphys.h"
82
83	REAL,INTENT(IN) :: time_0 ! not used
84
85	c dynamical variables:
86	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
87	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
88	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
89	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
90	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
91	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
92	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
93
94	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
95	REAL pks(ip1jmp1) ! exner at the surface
96	REAL pk(ip1jmp1,llm) ! exner at mid-layer
97	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
98	REAL phi(ip1jmp1,llm) ! geopotential
99	REAL w(ip1jmp1,llm) ! vertical velocity
100	REAL kpd(ip1jmp1) ! TB15 = exp (-z/H)
101
102	! ADAPTATION GCM POUR CP(T)
103	REAL temp(ip1jmp1,llm) ! temperature
104	REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk
105
106	! real zqmin,zqmax
107
108
109	! ED18 nogcm
110	REAL tau_ps
111	REAL tau_co2
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 pqco2mean ! moyenne globale ps*qco
202	REAL p0 ! pression de reference
203	REAL p00d ! globalaverage(kpd)
204	REAL qmean_co2,qmean_co2_vert ! mass mean mixing ratio vap co2
205	REAL pqco2(ip1jmp1) ! average co2 mass index : ps*q_co2
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 :: co2mass(iip1,jjp1)
244	real :: co2ice_ij(iip1,jjp1)
245	integer,save :: igcm_co2=0 ! index of CO2 tracer (if any)
246	integer :: i,j,ig
247	integer, parameter :: ngrid = 2+(jjm-1)*iim
248	! local mass
249	real :: mq(ip1jmp1,llm)
250
251	if (nday>=0) then
252	itaufin = nday*day_step
253	else
254	! to run a given (-nday) number of dynamical steps
255	itaufin = -nday
256	endif
257	if (less1day) then
258	c MODIF VENUS: to run less than one day:
259	itaufin = int(fractday*day_step)
260	endif
261	if (ndynstep.gt.0) then
262	! running a given number of dynamical steps
263	itaufin=ndynstep
264	endif
265	itaufinp1 = itaufin +1
266
267
268
269
270	c INITIALISATIONS
271	dudis(:,:) =0.
272	dvdis(:,:) =0.
273	dtetadis(:,:)=0.
274	dutop(:,:) =0.
275	c dvtop(:,:) =0.
276	c dtetatop(:,:)=0.
277	c dqtop(:,:,:) =0.
278	c dptop(:) =0.
279	dufi(:,:) =0.
280	dvfi(:,:) =0.
281	dtetafi(:,:)=0.
282	dqfi(:,:,:) =0.
283	dpfi(:) =0.
284
285	itau = 0
286	physics=.true.
287	if (iflag_phys==0.or.iflag_phys==2) physics=.false.
288
289	! ED18 nogcm
290	! firstcall_globaverage2d = .true.
291
292	c iday = day_ini+itau/day_step
293	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
294	c IF(time.GT.1.) THEN
295	c time = time-1.
296	c iday = iday+1
297	c ENDIF
298
299
300	c-----------------------------------------------------------------------
301	c On initialise la pression et la fonction d'Exner :
302	c --------------------------------------------------
303
304	dq(:,:,:)=0.
305	CALL pression ( ip1jmp1, ap, bp, ps, p )
306	if (pressure_exner) then
307	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
308	else
309	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
310	endif
311
312	c------------------
313	c TEST PK MONOTONE
314	c------------------
315	write(,) "Test PK"
316	do ij=1,ip1jmp1
317	do l=2,llm
318	! PRINT*,'pk(ij,l) = ',pk(ij,l)
319	! PRINT*,'pk(ij,l-1) = ',pk(ij,l-1)
320	if(pk(ij,l).gt.pk(ij,l-1)) then
321	c write(,) ij,l,pk(ij,l)
322	abort_message = 'PK non strictement decroissante'
323	call abort_gcm(modname,abort_message,1)
324	c write(,) "ATTENTION, Test PK deconnecte..."
325	endif
326	enddo
327	enddo
328	write(,) "Fin Test PK"
329	c stop
330	c------------------
331	c Preparing mixing of pressure and tracers in nogcm
332	c kpd=exp(-z/H) Needed to define a reference pressure :
333	c P0=pmean/globalaverage(kpd)
334	c thus P(i) = P0exp(-z(i)/H) = P0kpd(i)
335	c It is checked that globalaverage2d(Pi)=pmean
336	DO ij=1,ip1jmp1
337	kpd(ij) = exp(-phis(ij)/(r*200.))
338	ENDDO
339	p00d=globaverage2d(kpd) ! mean pres at ref level
340	tau_ps = 1. ! constante de rappel for pressure (s)
341	tau_co2 = 1.E5 !E5 ! constante de rappel for mix ratio qco2 (s)
342	tau_teta = 1.E7 !constante de rappel for potentiel temperature
343
344	! ED18 TEST
345	! PRINT*,'igcm_co2 = ',igcm_co2
346	! Locate tracer "co2" and set igcm_co2:
347	do iq=1,nqtot
348	if (tname(iq)=="co2") then
349	igcm_co2=iq
350	exit
351	endif
352	enddo
353
354	c-----------------------------------------------------------------------
355	c Debut de l'integration temporelle:
356	c ----------------------------------
357
358	1 CONTINUE ! Matsuno Forward step begins here
359
360	c date: (NB: date remains unchanged for Backward step)
361	c -----
362
363	jD_cur = jD_ref + day_ini - day_ref + &
364	& (itau+1)/day_step
365	jH_cur = jH_ref + start_time + &
366	& mod(itau+1,day_step)/float(day_step)
367	jD_cur = jD_cur + int(jH_cur)
368	jH_cur = jH_cur - int(jH_cur)
369
370	c
371
372	! Save fields obtained at previous time step as '...m1'
373	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
374	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
375	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
376	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
377	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
378
379	forward = .TRUE.
380	leapf = .FALSE.
381	dt = dtvr
382
383	2 CONTINUE ! Matsuno backward or leapfrog step begins here
384
385	c-----------------------------------------------------------------------
386
387	c date: (NB: only leapfrog step requires recomputing date)
388	c -----
389
390	IF (leapf) THEN
391	jD_cur = jD_ref + day_ini - day_ref +
392	& (itau+1)/day_step
393	jH_cur = jH_ref + start_time +
394	& mod(itau+1,day_step)/float(day_step)
395	jD_cur = jD_cur + int(jH_cur)
396	jH_cur = jH_cur - int(jH_cur)
397	ENDIF
398
399
400	c gestion des appels de la physique et des dissipations:
401	c ------------------------------------------------------
402	c
403	c ... P.Le Van ( 6/02/95 ) ....
404
405	! ED18: suppression des mentions de la variable apdiss dans le cas
406	! 'nogcm'
407
408	apphys = .FALSE.
409	statcl = .FALSE.
410	conser = .FALSE.
411
412
413	IF( purmats ) THEN
414	! Purely Matsuno time stepping
415	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
416
417
418	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
419	s .and. physics ) apphys = .TRUE.
420	ELSE
421	! Leapfrog/Matsuno time stepping
422	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
423
424
425	IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
426	END IF
427
428	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
429	! supress dissipation step
430
431
432
433	c-----------------------------------------------------------------------
434	c calcul des tendances dynamiques:
435	c --------------------------------
436
437	! ED18: suppression de l'onglet pour le cas nogcm
438
439
440	dv(:,:) = 0.
441	du(:,:) = 0.
442	dteta(:,:) = 0.
443	dq(:,:,:) = 0.
444
445
446
447	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
448	c
449	c-----------------------------------------------------------------------
450	c calcul des tendances physiques:
451	c -------------------------------
452	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
453	c
454	IF( purmats ) THEN
455	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
456	ELSE
457	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
458	ENDIF
459	c
460	c
461	IF( apphys ) THEN
462	c
463	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
464	c
465
466	CALL pression ( ip1jmp1, ap, bp, ps, p )
467	if (pressure_exner) then
468	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
469	else
470	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
471	endif
472
473	! Compute geopotential (physics might need it)
474	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
475
476	jD_cur = jD_ref + day_ini - day_ref + &
477	& (itau+1)/day_step
478
479	IF ((planet_type .eq."generic").or.
480	& (planet_type .eq."mars")) THEN
481	! AS: we make jD_cur to be pday
482	jD_cur = int(day_ini + itau/day_step)
483	ENDIF
484
485	! print*,'itau =',itau
486	! print*,'day_step =',day_step
487	! print*,'itau/day_step =',itau/day_step
488
489
490	jH_cur = jH_ref + start_time + &
491	& mod(itau+1,day_step)/float(day_step)
492	IF ((planet_type .eq."generic").or.
493	& (planet_type .eq."mars")) THEN
494	jH_cur = jH_ref + start_time + &
495	& mod(itau,day_step)/float(day_step)
496	ENDIF
497	jD_cur = jD_cur + int(jH_cur)
498	jH_cur = jH_cur - int(jH_cur)
499
500
501
502	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
503	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
504	! write(lunout,*)'current date = ',an, mois, jour, secondes
505
506	c rajout debug
507	c lafin = .true.
508
509
510	c Interface avec les routines de phylmd (phymars ... )
511	c -----------------------------------------------------
512
513	c+jld
514
515	c Diagnostique de conservation de l'Energie : initialisation
516	IF (ip_ebil_dyn.ge.1 ) THEN
517	ztit='bil dyn'
518	! Ehouarn: be careful, diagedyn is Earth-specific!
519	IF (planet_type.eq."earth") THEN
520	CALL diagedyn(ztit,2,1,1,dtphys
521	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
522	ENDIF
523	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
524	c-jld
525	#ifdef CPP_IOIPSL
526	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
527	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
528	c IF (first) THEN
529	c first=.false.
530	c#include "ini_paramLMDZ_dyn.h"
531	c ENDIF
532	c
533	c#include "write_paramLMDZ_dyn.h"
534	c
535	#endif
536	! #endif of #ifdef CPP_IOIPSL
537
538	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
539	! print*,'---------LEAPFROG---------------'
540	! print*,''
541	! print*,'> AVANT CALFIS'
542	! print*,''
543	! print*,'teta(3049,:) = ',teta(3049,:)
544	! print*,''
545	! print*,'dteta(3049,:) = ',dteta(3049,:)
546	! print*,''
547	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
548	! print*,''
549
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_co2)
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_co2) = 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_co2 with physical tendancy
660	! IF (igcm_co2.ne.0) then
661	! DO l=1, llm
662	! DO ij=1,ip1jmp1
663	! q(ij,l,igcm_co2)=q(ij,l,igcm_co2)+
664	! & dqfi(ij,l,igcm_co2)*dtphys
665	! ENDDO
666	! ENDDO
667	! ENDIF
668
669	c Mixing CO2 vertically
670	c --------------------------
671	if (igcm_co2.ne.0) then
672	DO ij=1,ip1jmp1
673	qmean_co2_vert=0.
674	DO l=1, llm
675	qmean_co2_vert= qmean_co2_vert
676	& + q(ij,l,igcm_co2)*( p(ij,l) - p(ij,l+1))
677	END DO
678	qmean_co2_vert= qmean_co2_vert/ps(ij)
679	DO l=1, llm
680	q(ij,l,igcm_co2)= qmean_co2_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 qco2
696	if (igcm_co2.ne.0) then
697
698	DO ij=1,ip1jmp1
699	pqco2(ij)= ps(ij) * q(ij,1,igcm_co2)
700	ENDDO
701	pqco2mean=globaverage2d(pqco2)
702
703	! Rappel newtonien vers qco2_mean
704	qmean_co2= pqco2mean / psmean
705
706	DO ij=1,ip1jmp1
707	q(ij,1,igcm_co2)=q(ij,1,igcm_co2)+
708	& (qmean_co2-q(ij,1,igcm_co2))*
709	& (1.-exp(-dtphys/tau_co2))
710	ENDDO
711
712	DO l=2, llm
713	DO ij=1,ip1jmp1
714	q(ij,l,igcm_co2)=q(ij,1,igcm_co2)
715	END DO
716	END DO
717
718	! TEMPORAIRE (ED)
719	! PRINT*,'psmean = ',psmean
720	! PRINT*,'qmean_co2 = ',qmean_co2
721	! PRINT*,'pqco2mean = ',pqco2mean
722	! PRINT*,'q(50,1,igcm_co2) = ',q(50,1,igcm_co2)
723	! PRINT*,'q(50,2,igcm_co2) = ',q(50,2,igcm_co2)
724	! PRINT*,'q(50,3,igcm_co2) = ',q(50,3,igcm_co2)
725
726	endif ! igcm_co2.ne.0
727
728
729	! ********************************************************
730
731	c Horizontal mixing of pressure
732	c ------------------------------
733	! Rappel newtonien vers psmean
734	psmean= globaverage2d(ps) ! mean pressure
735
736
737	c Mixing Temperature horizontally
738	c -------------------------------
739	! initialize variables that will be averaged
740	DO l=1,llm
741	DO ij=1,ip1jmp1
742	dp(ij,l) = p(ij,l) - p(ij,l+1)
743	tetadp(ij,l) = teta(ij,l)*dp(ij,l)
744	ENDDO
745	ENDDO
746
747	DO l=1,llm
748	tetadpmean = globaverage2d(tetadp(:,l))
749	dpmean = globaverage2d(dp(:,l))
750	tetamean = tetadpmean / dpmean
751	DO ij=1,ip1jmp1
752	teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) *
753	& (1 - exp(-dtphys/tau_teta))
754	ENDDO
755	ENDDO
756
757
758	ENDIF ! of IF( apphys )
759
760
761	c ********************************************************************
762	c ********************************************************************
763	c .... fin de l'integration dynamique et physique pour le pas itau ..
764	c ********************************************************************
765	c ********************************************************************
766
767
768	IF ( .NOT.purmats ) THEN
769	c ........................................................
770	c .............. schema matsuno + leapfrog ..............
771	c ........................................................
772
773	IF(forward. OR. leapf) THEN
774	itau= itau + 1
775	c iday= day_ini+itau/day_step
776	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
777	c IF(time.GT.1.) THEN
778	c time = time-1.
779	c iday = iday+1
780	c ENDIF
781	ENDIF
782
783
784	IF( itau. EQ. itaufinp1 ) then
785	if (flag_verif) then
786	write(79,*) 'ucov',ucov
787	write(80,*) 'vcov',vcov
788	write(81,*) 'teta',teta
789	write(82,*) 'ps',ps
790	write(83,*) 'q',q
791	WRITE(85,*) 'q1 = ',q(:,:,1)
792	WRITE(86,*) 'q3 = ',q(:,:,3)
793	endif
794
795	abort_message = 'Simulation finished'
796
797	call abort_gcm(modname,abort_message,0)
798	ENDIF
799	c-----------------------------------------------------------------------
800	c ecriture du fichier histoire moyenne:
801	c -------------------------------------
802
803	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
804	IF(itau.EQ.itaufin) THEN
805	iav=1
806	ELSE
807	iav=0
808	ENDIF
809
810	! ! Ehouarn: re-compute geopotential for outputs
811	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
812
813	IF (ok_dynzon) THEN
814	#ifdef CPP_IOIPSL
815	c les traceurs ne sont pas sortis, trop lourd.
816	c Peut changer eventuellement si besoin.
817	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
818	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
819	& du,dudis,dutop,dufi)
820	#endif
821	END IF
822	IF (ok_dyn_ave) THEN
823	#ifdef CPP_IOIPSL
824	CALL writedynav(itau,vcov,
825	& ucov,teta,pk,phi,q,masse,ps,phis)
826	#endif
827	ENDIF
828
829	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
830
831	if (ok_iso_verif) then
832	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
833	endif !if (ok_iso_verif) then
834
835	c-----------------------------------------------------------------------
836	c ecriture de la bande histoire:
837	c ------------------------------
838
839	IF( MOD(itau,iecri).EQ.0) THEN
840	! Ehouarn: output only during LF or Backward Matsuno
841	if (leapf.or.(.not.leapf.and.(.not.forward))) then
842	! ADAPTATION GCM POUR CP(T)
843	call tpot2t(ijp1llm,teta,temp,pk)
844	tsurpk = cpp*temp/pk
845	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
846	unat=0.
847	do l=1,llm
848	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
849	vnat(:,l)=vcov(:,l)/cv(:)
850	enddo
851	#ifdef CPP_IOIPSL
852	if (ok_dyn_ins) then
853	! write(lunout,*) "leapfrog: call writehist, itau=",itau
854	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
855	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
856	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
857	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
858	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
859	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
860	endif ! of if (ok_dyn_ins)
861	#endif
862	! For some Grads outputs of fields
863	if (output_grads_dyn) then
864	#include "write_grads_dyn.h"
865	endif
866	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
867	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
868
869	IF(itau.EQ.itaufin) THEN
870
871	if (planet_type=="mars") then
872	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
873	& vcov,ucov,teta,q,masse,ps)
874	else
875	CALL dynredem1("restart.nc",start_time,
876	& vcov,ucov,teta,q,masse,ps)
877	endif
878	CLOSE(99)
879	!!! Ehouarn: Why not stop here and now?
880	ENDIF ! of IF (itau.EQ.itaufin)
881
882	c-----------------------------------------------------------------------
883	c gestion de l'integration temporelle:
884	c ------------------------------------
885
886	IF( MOD(itau,iperiod).EQ.0 ) THEN
887	GO TO 1
888	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
889
890	IF( forward ) THEN
891	c fin du pas forward et debut du pas backward
892
893	forward = .FALSE.
894	leapf = .FALSE.
895	GO TO 2
896
897	ELSE
898	c fin du pas backward et debut du premier pas leapfrog
899
900	leapf = .TRUE.
901	dt = 2.*dtvr
902	GO TO 2
903	END IF ! of IF (forward)
904	ELSE
905
906	c ...... pas leapfrog .....
907
908	leapf = .TRUE.
909	dt = 2.*dtvr
910	GO TO 2
911	END IF ! of IF (MOD(itau,iperiod).EQ.0)
912	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
913
914	ELSE ! of IF (.not.purmats)
915
916	if (ok_iso_verif) then
917	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
918	endif !if (ok_iso_verif) then
919
920	c ........................................................
921	c .............. schema matsuno ...............
922	c ........................................................
923	IF( forward ) THEN
924
925	itau = itau + 1
926	c iday = day_ini+itau/day_step
927	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
928	c
929	c IF(time.GT.1.) THEN
930	c time = time-1.
931	c iday = iday+1
932	c ENDIF
933
934	forward = .FALSE.
935	IF( itau. EQ. itaufinp1 ) then
936	abort_message = 'Simulation finished'
937	call abort_gcm(modname,abort_message,0)
938	ENDIF
939	GO TO 2
940
941	ELSE ! of IF(forward) i.e. backward step
942
943	if (ok_iso_verif) then
944	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
945	endif !if (ok_iso_verif) then
946
947	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
948	IF(itau.EQ.itaufin) THEN
949	iav=1
950	ELSE
951	iav=0
952	ENDIF
953
954	! ! Ehouarn: re-compute geopotential for outputs
955	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
956
957	IF (ok_dynzon) THEN
958	#ifdef CPP_IOIPSL
959	c les traceurs ne sont pas sortis, trop lourd.
960	c Peut changer eventuellement si besoin.
961	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
962	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
963	& du,dudis,dutop,dufi)
964	#endif
965	ENDIF
966	IF (ok_dyn_ave) THEN
967	#ifdef CPP_IOIPSL
968	CALL writedynav(itau,vcov,
969	& ucov,teta,pk,phi,q,masse,ps,phis)
970	#endif
971	ENDIF
972
973	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
974
975	IF(MOD(itau,iecri ).EQ.0) THEN
976	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
977	! ADAPTATION GCM POUR CP(T)
978	call tpot2t(ijp1llm,teta,temp,pk)
979	tsurpk = cpp*temp/pk
980	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
981	unat=0.
982	do l=1,llm
983	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
984	vnat(:,l)=vcov(:,l)/cv(:)
985	enddo
986	#ifdef CPP_IOIPSL
987	if (ok_dyn_ins) then
988	! write(lunout,*) "leapfrog: call writehist (b)",
989	! & itau,iecri
990	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
991	endif ! of if (ok_dyn_ins)
992	#endif
993	! For some Grads outputs
994	if (output_grads_dyn) then
995	#include "write_grads_dyn.h"
996	endif
997
998	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
999
1000	IF(itau.EQ.itaufin) THEN
1001	if (planet_type=="mars") then
1002	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
1003	& vcov,ucov,teta,q,masse,ps)
1004	else
1005	CALL dynredem1("restart.nc",start_time,
1006	& vcov,ucov,teta,q,masse,ps)
1007	endif
1008	ENDIF ! of IF(itau.EQ.itaufin)
1009
1010	forward = .TRUE.
1011	GO TO 1
1012
1013	ENDIF ! of IF (forward)
1014
1015	END IF ! of IF(.not.purmats)
1016
1017	STOP
1018	END
1019
1020	! ************************************************************
1021	! globalaverage VERSION PLuto
1022	! ************************************************************
1023	FUNCTION globaverage2d(var)
1024	! ************************************************************
1025	IMPLICIT NONE
1026	#include "dimensions.h"
1027	#include "paramet.h"
1028	#include "comgeom2.h"
1029	!ip1jmp1 called in paramet.h = iip1 x jjp1
1030	REAL var(iip1,jjp1), globaverage2d
1031	integer i,j
1032	REAL airetot
1033	save airetot
1034	logical firstcall
1035	data firstcall/.true./
1036	save firstcall
1037
1038	if (firstcall) then
1039	firstcall=.false.
1040	airetot =0.
1041	DO j=2,jjp1-1 ! lat
1042	DO i = 1, iim ! lon
1043	airetot = airetot + aire(i,j)
1044	END DO
1045	END DO
1046	DO i=1,iim
1047	airetot=airetot + aire(i,1)+aire(i,jjp1)
1048	ENDDO
1049	end if
1050
1051	globaverage2d = 0.
1052	DO j=2,jjp1-1
1053	DO i = 1, iim
1054	globaverage2d = globaverage2d + var(i,j)*aire(i,j)
1055	END DO
1056	END DO
1057
1058	DO i=1,iim
1059	globaverage2d = globaverage2d + var(i,1)*aire(i,1)
1060	globaverage2d = globaverage2d + var(i,jjp1)*aire(i,jjp1)
1061	ENDDO
1062
1063	globaverage2d = globaverage2d/airetot
1064	return
1065	end
1066

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

Download in other formats:

Original Format