Context Navigation

advtrac_loc.F @ 2337

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

fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
comment out many invasive debug writes.

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

File size: 13.0 KB

Rev	Line
[1632]	1	!
[1673]	2	! $Id$
[1632]	3	!
	4	c
	5	c
	6	#define DEBUG_IO
	7	#undef DEBUG_IO
	8	SUBROUTINE advtrac_loc(pbarug,pbarvg ,wg,
	9	* p, massem,q,teta,
	10	* pk )
	11
	12	c Auteur : F. Hourdin
	13	c
	14	c Modif. P. Le Van (20/12/97)
	15	c F. Codron (10/99)
	16	c D. Le Croller (07/2001)
	17	c M.A Filiberti (04/2002)
	18	c
[1823]	19	USE parallel_lmdz
[1632]	20	USE Write_Field_loc
	21	USE Write_Field
	22	USE Bands
	23	USE mod_hallo
	24	USE Vampir
	25	USE times
[2270]	26	USE infotrac, ONLY: nqtot, iadv, ok_iso_verif
[1987]	27	USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
	28	USE advtrac_mod, ONLY: finmasse
[1632]	29	IMPLICIT NONE
	30	c
	31	#include "dimensions.h"
	32	#include "paramet.h"
	33	#include "comconst.h"
	34	#include "comvert.h"
	35	#include "comdissip.h"
	36	#include "comgeom2.h"
	37	#include "logic.h"
	38	#include "temps.h"
	39	#include "ener.h"
	40	#include "description.h"
	41
	42	c-------------------------------------------------------------------
	43	c Arguments
	44	c-------------------------------------------------------------------
	45	c Ajout PPM
	46	c--------------------------------------------------------
	47	REAL massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
	48	c--------------------------------------------------------
	49	INTEGER iapptrac
	50	REAL pbarug(ijb_u:ije_u,llm),pbarvg(ijb_v:ije_v,llm)
	51	REAL wg(ijb_u:ije_u,llm)
	52	REAL q(ijb_u:ije_u,llm,nqtot),massem(ijb_u:ije_u,llm)
	53	REAL p( ijb_u:ije_u,llmp1 ),teta(ijb_u:ije_u,llm)
	54	REAL pk(ijb_u:ije_u,llm)
	55
	56	c-------------------------------------------------------------
	57	c Variables locales
	58	c-------------------------------------------------------------
	59
	60	REAL zdp(ijb_u:ije_u)
	61	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
	62	INTEGER,SAVE :: iadvtr=0
	63	c$OMP THREADPRIVATE(iadvtr)
	64	INTEGER ij,l,iq,iiq
	65	REAL zdpmin, zdpmax
	66	c----------------------------------------------------------
	67	c Rajouts pour PPM
	68	c----------------------------------------------------------
	69	INTEGER indice,n
	70	REAL dtbon ! Pas de temps adaptatif pour que CFL<1
	71	REAL CFLmaxz,aaa,bbb ! CFL maximum
	72	REAL psppm(iim,jjb_u:jje_u) ! pression au sol
	73	REAL unatppm(iim,jjb_u:jje_u,llm),vnatppm(iim,jjb_u:jje_u,llm)
	74	REAL qppm(iim*jjnb_u,llm,nqtot)
	75	REAL fluxwppm(iim,jjb_u:jje_u,llm)
	76	REAL apppm(llmp1), bpppm(llmp1)
	77	LOGICAL dum,fill
	78	DATA fill/.true./
	79	DATA dum/.true./
	80	integer ijb,ije,ijbu,ijbv,ijeu,ijev,j
[1848]	81	type(Request),SAVE :: testRequest
	82	!$OMP THREADPRIVATE(testRequest)
[1632]	83
[2270]	84	c test sur l''eventuelle creation de valeurs negatives de la masse
[1632]	85	ijb=ij_begin
	86	ije=ij_end
	87	if (pole_nord) ijb=ij_begin+iip1
	88	if (pole_sud) ije=ij_end-iip1
	89
	90	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	91	DO l=1,llm-1
	92	DO ij = ijb+1,ije
	93	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)
	94	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
	95	s + wg(ij,l+1) - wg(ij,l)
	96	ENDDO
	97
	98	c CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
	99	c ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
	100
	101	do ij=ijb,ije-iip1+1,iip1
	102	zdp(ij)=zdp(ij+iip1-1)
	103	enddo
	104
	105	DO ij = ijb,ije
	106	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
	107	ENDDO
	108
	109
	110	c CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
	111	c ym ---> eventuellement a revoir
	112	CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
	113
	114	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
	115	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,
	116	s ' MAX:', zdpmax
	117	ENDIF
	118
	119	ENDDO
	120	c$OMP END DO NOWAIT
	121
	122	c-------------------------------------------------------------------
	123	c Advection proprement dite (Modification Le Croller (07/2001)
	124	c-------------------------------------------------------------------
	125
	126	c----------------------------------------------------
	127	c Calcul des moyennes basées sur la masse
	128	c----------------------------------------------------
	129
	130	cym ----> Normalement, inutile pour les schémas classiques
	131	cym ----> Revérifier lors de la parallélisation des autres schemas
	132
	133	cym call massbar_p(massem,massebx,masseby)
	134
	135	#ifdef DEBUG_IO
	136	CALL WriteField_u('massem',massem)
	137	CALL WriteField_u('wg',wg)
	138	CALL WriteField_u('pbarug',pbarug)
	139	CALL WriteField_v('pbarvg',pbarvg)
	140	CALL WriteField_u('p_tmp',p)
	141	CALL WriteField_u('pk_tmp',pk)
	142	CALL WriteField_u('teta_tmp',teta)
	143	do j=1,nqtot
	144	call WriteField_u('q_adv'//trim(int2str(j)),
	145	. q(:,:,j))
	146	enddo
	147	#endif
	148
	149	!
	150	! call Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
	151	!
	152	! call SendRequest(TestRequest)
	153	!c$OMP BARRIER
	154	! call WaitRequest(TestRequest)
	155	c$OMP BARRIER
	156
[2286]	157	!write(,) 'advtrac 157: appel de vlspltgen_loc'
[1632]	158	call vlspltgen_loc( q,iadv, 2., massem, wg ,
	159	* pbarug,pbarvg,dtvr,p,
	160	* pk,teta )
	161
[2286]	162	!write(,) 'advtrac 162: apres appel vlspltgen_loc'
[2270]	163	if (ok_iso_verif) then
	164	call check_isotopes(q,ijb_u,ije_u,'advtrac 162')
	165	endif !if (ok_iso_verif) then
	166
[1632]	167	#ifdef DEBUG_IO
	168	do j=1,nqtot
	169	call WriteField_u('q_adv'//trim(int2str(j)),
	170	. q(:,:,j))
	171	enddo
	172	#endif
	173
	174	GOTO 1234
	175	c-----------------------------------------------------------
	176	c Appel des sous programmes d'advection
	177	c-----------------------------------------------------------
	178	do iq=1,nqtot
	179	c call clock(t_initial)
	180	if(iadv(iq) == 0) cycle
	181	c ----------------------------------------------------------------
	182	c Schema de Van Leer I MUSCL
	183	c ----------------------------------------------------------------
	184	if(iadv(iq).eq.10) THEN
	185
	186	!LF call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
	187
	188	c ----------------------------------------------------------------
	189	c Schema "pseudo amont" + test sur humidite specifique
	190	C pour la vapeur d'eau. F. Codron
	191	c ----------------------------------------------------------------
	192	else if(iadv(iq).eq.14) then
	193	c
	194	cym stop 'advtrac : appel à vlspltqs :schema non parallelise'
	195	!LF CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
	196	!LF * pbarug,pbarvg,dtvr,p,pk,teta )
	197	c ----------------------------------------------------------------
	198	c Schema de Frederic Hourdin
	199	c ----------------------------------------------------------------
	200	else if(iadv(iq).eq.12) then
	201	stop 'advtrac : schema non parallelise'
	202	c Pas de temps adaptatif
	203	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	204	if (n.GT.1) then
	205	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	206	s dtvr,'n=',n
	207	endif
	208	do indice=1,n
	209	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
	210	end do
	211	else if(iadv(iq).eq.13) then
	212	stop 'advtrac : schema non parallelise'
	213	c Pas de temps adaptatif
	214	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	215	if (n.GT.1) then
	216	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	217	s dtvr,'n=',n
	218	endif
	219	do indice=1,n
	220	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
	221	end do
	222	c ----------------------------------------------------------------
	223	c Schema de pente SLOPES
	224	c ----------------------------------------------------------------
	225	else if (iadv(iq).eq.20) then
	226	stop 'advtrac : schema non parallelise'
	227
	228	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
	229
	230	c ----------------------------------------------------------------
	231	c Schema de Prather
	232	c ----------------------------------------------------------------
	233	else if (iadv(iq).eq.30) then
	234	stop 'advtrac : schema non parallelise'
	235	c Pas de temps adaptatif
	236	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	237	if (n.GT.1) then
	238	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	239	s dtvr,'n=',n
	240	endif
	241	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
	242	s n,dtbon)
	243	c ----------------------------------------------------------------
	244	c Schemas PPM Lin et Rood
	245	c ----------------------------------------------------------------
	246	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
	247	s iadv(iq).LE.18)) then
	248
	249	stop 'advtrac : schema non parallelise'
	250
	251	c Test sur le flux horizontal
	252	c Pas de temps adaptatif
	253	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
	254	if (n.GT.1) then
	255	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
	256	s dtvr,'n=',n
	257	endif
	258	c Test sur le flux vertical
	259	CFLmaxz=0.
	260	do l=2,llm
	261	do ij=iip2,ip1jm
	262	aaa=wg(ij,l)*dtvr/massem(ij,l)
	263	CFLmaxz=max(CFLmaxz,aaa)
	264	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
	265	CFLmaxz=max(CFLmaxz,bbb)
	266	enddo
	267	enddo
	268	if (CFLmaxz.GE.1) then
	269	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
	270	endif
	271
	272	c-----------------------------------------------------------
	273	c Ss-prg interface LMDZ.4->PPM3d
	274	c-----------------------------------------------------------
	275
	276	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
	277	s apppm,bpppm,massebx,masseby,pbarug,pbarvg,
	278	s unatppm,vnatppm,psppm)
	279
	280	do indice=1,n
	281	c---------------------------------------------------------------------
	282	c VL (version PPM) horiz. et PPM vert.
	283	c---------------------------------------------------------------------
	284	if (iadv(iq).eq.11) then
	285	c Ss-prg PPM3d de Lin
	286	call ppm3d(1,qppm(1,1,iq),
	287	s psppm,psppm,
	288	s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
	289	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	290	s fill,dum,220.)
	291
	292	c----------------------------------------------------------------------
	293	c Monotonic PPM
	294	c----------------------------------------------------------------------
	295	else if (iadv(iq).eq.16) then
	296	c Ss-prg PPM3d de Lin
	297	call ppm3d(1,qppm(1,1,iq),
	298	s psppm,psppm,
	299	s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
	300	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	301	s fill,dum,220.)
	302	c---------------------------------------------------------------------
	303
	304	c---------------------------------------------------------------------
	305	c Semi Monotonic PPM
	306	c---------------------------------------------------------------------
	307	else if (iadv(iq).eq.17) then
	308	c Ss-prg PPM3d de Lin
	309	call ppm3d(1,qppm(1,1,iq),
	310	s psppm,psppm,
	311	s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
	312	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	313	s fill,dum,220.)
	314	c---------------------------------------------------------------------
	315
	316	c---------------------------------------------------------------------
	317	c Positive Definite PPM
	318	c---------------------------------------------------------------------
	319	else if (iadv(iq).eq.18) then
	320	c Ss-prg PPM3d de Lin
	321	call ppm3d(1,qppm(1,1,iq),
	322	s psppm,psppm,
	323	s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
	324	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
	325	s fill,dum,220.)
	326	c---------------------------------------------------------------------
	327	endif
	328	enddo
	329	c-----------------------------------------------------------------
	330	c Ss-prg interface PPM3d-LMDZ.4
	331	c-----------------------------------------------------------------
	332	call interpost(q(1,1,iq),qppm(1,1,iq))
	333	endif
	334	c----------------------------------------------------------------------
	335
	336	c-----------------------------------------------------------------
	337	c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
	338	c et Nord j=1
	339	c-----------------------------------------------------------------
	340
	341	c call traceurpole(q(1,1,iq),massem)
	342
	343	c calcul du temps cpu pour un schema donne
	344
	345
	346	end DO
	347
	348	1234 CONTINUE
	349	c$OMP BARRIER
	350
[1673]	351	if (planet_type=="earth") then
	352
[1632]	353	ijb=ij_begin
	354	ije=ij_end
	355
	356	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	357	DO l = 1, llm
	358	DO ij = ijb, ije
	359	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
	360	ENDDO
	361	ENDDO
	362	c$OMP END DO
	363
[2270]	364	! CRisi: on passe nqtot et non nq
	365	CALL qminimum_loc( q, nqtot, finmasse )
[1632]	366
[1673]	367	endif ! of if (planet_type=="earth")
[1632]	368
	369	RETURN
	370	END
	371

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

Download in other formats:

Original Format