Context Navigation

advtrac_p.F @ 5423

Last change on this file since 5423 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 16.5 KB

Rev	Line
[630]	1	!
[1299]	2	! $Id: advtrac_p.F 1299 2010-01-20 14:27:21Z evignon $
[630]	3	!
	4	c
	5	c
	6	SUBROUTINE advtrac_p(pbaru,pbarv ,
	7	* p, masse,q,iapptrac,teta,
	8	* flxw,
[764]	9	* pk )
[630]	10
	11	c Auteur : F. Hourdin
	12	c
	13	c Modif. P. Le Van (20/12/97)
	14	c F. Codron (10/99)
	15	c D. Le Croller (07/2001)
	16	c M.A Filiberti (04/2002)
	17	c
	18	USE parallel
[764]	19	USE Write_Field_p
[630]	20	USE Bands
	21	USE mod_hallo
	22	USE Vampir
	23	USE times
[1146]	24	USE infotrac
[1299]	25	USE control_mod
[630]	26	IMPLICIT NONE
	27	c
	28	#include "dimensions.h"
	29	#include "paramet.h"
	30	#include "comconst.h"
	31	#include "comvert.h"
	32	#include "comdissip.h"
	33	#include "comgeom2.h"
	34	#include "logic.h"
	35	#include "temps.h"
	36	#include "ener.h"
	37	#include "description.h"
	38
	39	c-------------------------------------------------------------------
	40	c Arguments
	41	c-------------------------------------------------------------------
	42	c Ajout PPM
	43	c--------------------------------------------------------
	44	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
	45	c--------------------------------------------------------
	46	INTEGER iapptrac
	47	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
[1146]	48	REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
[630]	49	REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
	50	REAL pk(ip1jmp1,llm)
[960]	51	REAL :: flxw(ip1jmp1,llm)
[630]	52
	53	c-------------------------------------------------------------
	54	c Variables locales
	55	c-------------------------------------------------------------
	56
	57	REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
	58	REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
[764]	59	REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
[630]	60	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
	61	INTEGER iadvtr
	62	INTEGER ij,l,iq,iiq
	63	REAL zdpmin, zdpmax
	64	SAVE iadvtr, massem, pbaruc, pbarvc
	65	DATA iadvtr/0/
[764]	66	c$OMP THREADPRIVATE(iadvtr)
[630]	67	c----------------------------------------------------------
	68	c Rajouts pour PPM
	69	c----------------------------------------------------------
	70	INTEGER indice,n
	71	REAL dtbon ! Pas de temps adaptatif pour que CFL<1
	72	REAL CFLmaxz,aaa,bbb ! CFL maximum
	73	REAL psppm(iim,jjp1) ! pression au sol
	74	REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
[1146]	75	REAL qppm(iim*jjp1,llm,nqtot)
[630]	76	REAL fluxwppm(iim,jjp1,llm)
	77	REAL apppm(llmp1), bpppm(llmp1)
	78	LOGICAL dum,fill
	79	DATA fill/.true./
	80	DATA dum/.true./
[985]	81	REAL,SAVE :: finmasse(ip1jmp1,llm)
[630]	82	integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
	83	type(Request) :: Request_vanleer
[764]	84	REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
	85	REAL,SAVE :: teta_tmp(ip1jmp1,llm)
	86	REAL,SAVE :: pk_tmp(ip1jmp1,llm)
[1146]	87
[630]	88	ijb_u=ij_begin
	89	ije_u=ij_end
	90
	91	ijb_v=ij_begin-iip1
	92	ije_v=ij_end
	93	if (pole_nord) ijb_v=ij_begin
	94	if (pole_sud) ije_v=ij_end-iip1
	95
	96	IF(iadvtr.EQ.0) THEN
	97	c CALL initial0(ijp1llm,pbaruc)
	98	c CALL initial0(ijmllm,pbarvc)
[764]	99	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	100	DO l=1,llm
	101	pbaruc(ijb_u:ije_u,l)=0.
	102	pbarvc(ijb_v:ije_v,l)=0.
	103	ENDDO
	104	c$OMP END DO NOWAIT
[630]	105	ENDIF
	106
	107	c accumulation des flux de masse horizontaux
[764]	108	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	109	DO l=1,llm
	110	DO ij = ijb_u,ije_u
	111	pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
	112	ENDDO
	113	DO ij = ijb_v,ije_v
	114	pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
	115	ENDDO
	116	ENDDO
[764]	117	c$OMP END DO NOWAIT
[630]	118
	119	c selection de la masse instantannee des mailles avant le transport.
	120	IF(iadvtr.EQ.0) THEN
	121
	122	c CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
	123	ijb=ij_begin
	124	ije=ij_end
	125
[764]	126	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	127	DO l=1,llm
	128	massem(ijb:ije,l)=masse(ijb:ije,l)
	129	ENDDO
	130	c$OMP END DO NOWAIT
	131
[630]	132	ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
	133	c
	134	ENDIF
	135
	136	iadvtr = iadvtr+1
[764]	137
	138	c$OMP MASTER
[630]	139	iapptrac = iadvtr
[764]	140	c$OMP END MASTER
[630]	141
[764]	142	c Test pour savoir si on advecte a ce pas de temps
[630]	143
	144	IF ( iadvtr.EQ.iapp_tracvl ) THEN
[764]	145	c$OMP MASTER
[630]	146	call suspend_timer(timer_caldyn)
[764]	147	c$OMP END MASTER
[630]	148
	149	ijb=ij_begin
	150	ije=ij_end
	151
	152
	153	cc .. Modif P.Le Van ( 20/12/97 ) ....
	154	cc
	155
	156	c traitement des flux de masse avant advection.
	157	c 1. calcul de w
	158	c 2. groupement des mailles pres du pole.
	159
	160	CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
[985]	161
[764]	162	c$OMP BARRIER
[630]	163
[985]	164	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	165	DO l=1,llmp1
	166	p_tmp(ijb:ije,l)=p(ijb:ije,l)
	167	ENDDO
	168	c$OMP END DO NOWAIT
	169
	170	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	171	DO l=1,llm
	172	pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
	173	teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
	174	ENDDO
	175	c$OMP END DO NOWAIT
	176
	177	c$OMP MASTER
[630]	178	call VTb(VTHallo)
[985]	179	c$OMP END MASTER
	180
[630]	181	call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm,
	182	* jj_Nb_vanleer,0,0,Request_vanleer)
	183	call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm,
	184	* jj_Nb_vanleer,1,0,Request_vanleer)
	185	call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm,
	186	* jj_Nb_vanleer,0,0,Request_vanleer)
	187	call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm,
	188	* jj_Nb_vanleer,0,0,Request_vanleer)
	189	call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm,
	190	* jj_Nb_vanleer,1,1,Request_vanleer)
	191	call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1,
	192	* jj_Nb_vanleer,1,1,Request_vanleer)
	193	call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,
	194	* jj_Nb_vanleer,1,1,Request_vanleer)
[1146]	195	do j=1,nqtot
[630]	196	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
	197	* jj_nb_vanleer,0,0,Request_vanleer)
	198	enddo
[985]	199
[630]	200	call SendRequest(Request_vanleer)
[985]	201	c$OMP BARRIER
[630]	202	call WaitRequest(Request_vanleer)
	203
[985]	204
	205	c$OMP BARRIER
	206	c$OMP MASTER
	207	call SetDistrib(jj_nb_vanleer)
[630]	208	call VTe(VTHallo)
	209	call VTb(VTadvection)
	210	call start_timer(timer_vanleer)
[764]	211	c$OMP END MASTER
	212	c$OMP BARRIER
[630]	213
	214	! ... Flux de masse diaganostiques traceurs
	215	ijb=ij_begin
	216	ije=ij_end
[1299]	217	flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
[630]	218
	219	c test sur l'eventuelle creation de valeurs negatives de la masse
	220	ijb=ij_begin
	221	ije=ij_end
	222	if (pole_nord) ijb=ij_begin+iip1
	223	if (pole_sud) ije=ij_end-iip1
	224
[764]	225	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	226	DO l=1,llm-1
	227	DO ij = ijb+1,ije
	228	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)
	229	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
	230	s + wg(ij,l+1) - wg(ij,l)
	231	ENDDO
	232
	233	c CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
	234	c ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
	235
	236	do ij=ijb,ije-iip1+1,iip1
	237	zdp(ij)=zdp(ij+iip1-1)
	238	enddo
	239
	240	DO ij = ijb,ije
	241	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
	242	ENDDO
	243
	244
	245	c CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
	246	c ym ---> eventuellement a revoir
	247	CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
	248
	249	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
	250	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,
	251	s ' MAX:', zdpmax
	252	ENDIF
	253
	254	ENDDO
[764]	255	c$OMP END DO NOWAIT
[630]	256
	257	c-------------------------------------------------------------------
	258	c Advection proprement dite (Modification Le Croller (07/2001)
	259	c-------------------------------------------------------------------
	260
	261	c----------------------------------------------------
	262	c Calcul des moyennes basées sur la masse
	263	c----------------------------------------------------
	264
	265	cym ----> Normalement, inutile pour les schémas classiques
	266	cym ----> Revérifier lors de la parallélisation des autres schemas
	267
[764]	268	cym call massbar_p(massem,massebx,masseby)
[630]	269
[764]	270	call vlspltgen_p( q,iadv, 2., massem, wg ,
	271	* pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
	272
	273
	274	GOTO 1234
[630]	275	c-----------------------------------------------------------
	276	c Appel des sous programmes d'advection
	277	c-----------------------------------------------------------
[1146]	278	do iq=1,nqtot
[630]	279	c call clock(t_initial)
	280	if(iadv(iq) == 0) cycle
	281	c ----------------------------------------------------------------
	282	c Schema de Van Leer I MUSCL
	283	c ----------------------------------------------------------------
	284	if(iadv(iq).eq.10) THEN
	285
	286	call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
	287
	288	c ----------------------------------------------------------------
	289	c Schema "pseudo amont" + test sur humidite specifique
	290	C pour la vapeur d'eau. F. Codron
	291	c ----------------------------------------------------------------
	292	else if(iadv(iq).eq.14) then
	293	c
	294	cym stop 'advtrac : appel à vlspltqs :schema non parallelise'
	295	CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
	296	* pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
	297	c ----------------------------------------------------------------
	298	c Schema de Frederic Hourdin
	299	c ----------------------------------------------------------------
	300	else if(iadv(iq).eq.12) then
	301	stop 'advtrac : schema non parallelise'
	302	c Pas de temps adaptatif
	303	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	304	if (n.GT.1) then
	305	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	306	s dtvr,'n=',n
	307	endif
	308	do indice=1,n
	309	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
	310	end do
	311	else if(iadv(iq).eq.13) then
	312	stop 'advtrac : schema non parallelise'
	313	c Pas de temps adaptatif
	314	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	315	if (n.GT.1) then
	316	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	317	s dtvr,'n=',n
	318	endif
	319	do indice=1,n
	320	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
	321	end do
	322	c ----------------------------------------------------------------
	323	c Schema de pente SLOPES
	324	c ----------------------------------------------------------------
	325	else if (iadv(iq).eq.20) then
	326	stop 'advtrac : schema non parallelise'
	327
	328	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
	329
	330	c ----------------------------------------------------------------
	331	c Schema de Prather
	332	c ----------------------------------------------------------------
	333	else if (iadv(iq).eq.30) then
	334	stop 'advtrac : schema non parallelise'
	335	c Pas de temps adaptatif
	336	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	337	if (n.GT.1) then
	338	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	339	s dtvr,'n=',n
	340	endif
	341	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
	342	s n,dtbon)
	343	c ----------------------------------------------------------------
	344	c Schemas PPM Lin et Rood
	345	c ----------------------------------------------------------------
	346	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
	347	s iadv(iq).LE.18)) then
	348
	349	stop 'advtrac : schema non parallelise'
	350
	351	c Test sur le flux horizontal
	352	c Pas de temps adaptatif
	353	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	354	if (n.GT.1) then
	355	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	356	s dtvr,'n=',n
	357	endif
	358	c Test sur le flux vertical
	359	CFLmaxz=0.
	360	do l=2,llm
	361	do ij=iip2,ip1jm
	362	aaa=wg(ij,l)*dtvr/massem(ij,l)
	363	CFLmaxz=max(CFLmaxz,aaa)
	364	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
	365	CFLmaxz=max(CFLmaxz,bbb)
	366	enddo
	367	enddo
	368	if (CFLmaxz.GE.1) then
	369	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
	370	endif
	371
	372	c-----------------------------------------------------------
	373	c Ss-prg interface LMDZ.4->PPM3d
	374	c-----------------------------------------------------------
	375
	376	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
	377	s apppm,bpppm,massebx,masseby,pbarug,pbarvg,
	378	s unatppm,vnatppm,psppm)
	379
	380	do indice=1,n
	381	c---------------------------------------------------------------------
	382	c VL (version PPM) horiz. et PPM vert.
	383	c---------------------------------------------------------------------
	384	if (iadv(iq).eq.11) then
	385	c Ss-prg PPM3d de Lin
	386	call ppm3d(1,qppm(1,1,iq),
	387	s psppm,psppm,
	388	s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
	389	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	390	s fill,dum,220.)
	391
	392	c----------------------------------------------------------------------
	393	c Monotonic PPM
	394	c----------------------------------------------------------------------
	395	else if (iadv(iq).eq.16) then
	396	c Ss-prg PPM3d de Lin
	397	call ppm3d(1,qppm(1,1,iq),
	398	s psppm,psppm,
	399	s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
	400	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	401	s fill,dum,220.)
	402	c---------------------------------------------------------------------
	403
	404	c---------------------------------------------------------------------
	405	c Semi Monotonic PPM
	406	c---------------------------------------------------------------------
	407	else if (iadv(iq).eq.17) then
	408	c Ss-prg PPM3d de Lin
	409	call ppm3d(1,qppm(1,1,iq),
	410	s psppm,psppm,
	411	s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
	412	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	413	s fill,dum,220.)
	414	c---------------------------------------------------------------------
	415
	416	c---------------------------------------------------------------------
	417	c Positive Definite PPM
	418	c---------------------------------------------------------------------
	419	else if (iadv(iq).eq.18) then
	420	c Ss-prg PPM3d de Lin
	421	call ppm3d(1,qppm(1,1,iq),
	422	s psppm,psppm,
	423	s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
	424	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	425	s fill,dum,220.)
	426	c---------------------------------------------------------------------
	427	endif
	428	enddo
	429	c-----------------------------------------------------------------
	430	c Ss-prg interface PPM3d-LMDZ.4
	431	c-----------------------------------------------------------------
	432	call interpost(q(1,1,iq),qppm(1,1,iq))
	433	endif
	434	c----------------------------------------------------------------------
	435
	436	c-----------------------------------------------------------------
	437	c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
	438	c et Nord j=1
	439	c-----------------------------------------------------------------
	440
	441	c call traceurpole(q(1,1,iq),massem)
	442
	443	c calcul du temps cpu pour un schema donne
	444
	445	c call clock(t_final)
	446	cym tps_cpu=t_final-t_initial
	447	cym cpuadv(iq)=cpuadv(iq)+tps_cpu
	448
	449	end DO
	450
[764]	451	1234 CONTINUE
	452	c$OMP BARRIER
[985]	453
[630]	454	ijb=ij_begin
	455	ije=ij_end
[985]	456
	457	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	458	DO l = 1, llm
	459	DO ij = ijb, ije
	460	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
	461	ENDDO
	462	ENDDO
[985]	463	c$OMP END DO
[630]	464
	465	CALL qminimum_p( q, 2, finmasse )
	466
	467	c------------------------------------------------------------------
	468	c on reinitialise a zero les flux de masse cumules
	469	c---------------------------------------------------
[764]	470	c iadvtr=0
[985]	471
	472	c$OMP MASTER
[630]	473	call VTe(VTadvection)
	474	call stop_timer(timer_vanleer)
[985]	475	call VTb(VThallo)
	476	c$OMP END MASTER
[630]	477
[1146]	478	do j=1,nqtot
[630]	479	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
	480	* jj_nb_caldyn,0,0,Request_vanleer)
	481	enddo
	482
[960]	483	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
	484	* jj_nb_caldyn,0,0,Request_vanleer)
	485
[630]	486	call SendRequest(Request_vanleer)
[985]	487	c$OMP BARRIER
[630]	488	call WaitRequest(Request_vanleer)
[985]	489
	490	c$OMP BARRIER
	491	c$OMP MASTER
	492	call SetDistrib(jj_nb_caldyn)
[630]	493	call VTe(VThallo)
	494	call resume_timer(timer_caldyn)
[764]	495	c$OMP END MASTER
	496	c$OMP BARRIER
	497	iadvtr=0
[630]	498	ENDIF ! if iadvtr.EQ.iapp_tracvl
	499
	500	RETURN
	501	END
	502

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

Download in other formats:

Original Format