Context Navigation

leapfrog_nogcm.F @ 3302

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

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

Download in other formats:

Original Format