Context Navigation

advtrac_p.F90 @ 2212

Last change on this file since 2212 was 1987, checked in by Ehouarn Millour, 10 years ago
Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency). Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid). EM
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.0 KB

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

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

Download in other formats:

Original Format