Context Navigation

advtrac_p.F @ 3900

Last change on this file since 3900 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

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

File size: 16.6 KB

Rev	Line
[630]	1	!
[1403]	2	! $Id: advtrac_p.F 1446 2010-10-22 09:27:25Z 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
[1403]	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
[1446]	134	ENDIF ! of IF(iadvtr.EQ.0)
[630]	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
[1403]	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
[1446]	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
[1446]	454	if (planet_type=="earth") then
[985]	455
[1446]	456	ijb=ij_begin
	457	ije=ij_end
	458
[985]	459	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1446]	460	DO l = 1, llm
[630]	461	DO ij = ijb, ije
	462	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
	463	ENDDO
[1446]	464	ENDDO
[985]	465	c$OMP END DO
[630]	466
[1446]	467	CALL qminimum_p( q, 2, finmasse )
[630]	468
	469	c------------------------------------------------------------------
	470	c on reinitialise a zero les flux de masse cumules
	471	c---------------------------------------------------
[764]	472	c iadvtr=0
[985]	473
	474	c$OMP MASTER
[630]	475	call VTe(VTadvection)
	476	call stop_timer(timer_vanleer)
[985]	477	call VTb(VThallo)
	478	c$OMP END MASTER
[630]	479
[1146]	480	do j=1,nqtot
[630]	481	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
	482	* jj_nb_caldyn,0,0,Request_vanleer)
	483	enddo
	484
[960]	485	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
	486	* jj_nb_caldyn,0,0,Request_vanleer)
	487
[630]	488	call SendRequest(Request_vanleer)
[985]	489	c$OMP BARRIER
[630]	490	call WaitRequest(Request_vanleer)
[985]	491
	492	c$OMP BARRIER
	493	c$OMP MASTER
	494	call SetDistrib(jj_nb_caldyn)
[630]	495	call VTe(VThallo)
	496	call resume_timer(timer_caldyn)
[764]	497	c$OMP END MASTER
	498	c$OMP BARRIER
	499	iadvtr=0
[1446]	500	endif ! of if (planet_type=="earth")
[630]	501	ENDIF ! if iadvtr.EQ.iapp_tracvl
	502
	503	RETURN
	504	END
	505

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

Download in other formats:

Original Format