Context Navigation

advtrac_p.F90 @ 2502

Last change on this file since 2502 was 2239, checked in by Ehouarn Millour, 10 years ago

Reorganizing physics/dynamics interface:

what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
"bibio" directory is now "misc" and should only contain autonomous utilities.
"cosp" is now a subdirectory of phylmd.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 17.2 KB

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

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

Download in other formats:

Original Format