Context Navigation

advtrac_loc.f90 @ 5748

Last change on this file since 5748 was 5748, checked in by dcugnet, 6 days ago

Use REAL(KIND=REAL32) and REAL(KIND=REAL64) Iinstead of REAL and DOUBLE PRECISION

to avoid ambiguity problems in generic procedure when reals are promoted to doubles.

generic "num2str" replaces "str2int", "str2real", "str2dble" and "str2bool" functions.

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: 12.2 KB

Rev	Line
[4056]	1
	2	SUBROUTINE advtrac_loc(pbarug, pbarvg, wg, p, massem, q, teta, pk)
	3	! Auteur : F. Hourdin
	4	!
	5	! Modif. P. Le Van (20/12/97)
	6	! F. Codron (10/99)
	7	! D. Le Croller (07/2001)
	8	! M.A Filiberti (04/2002)
	9	!
[5281]	10	USE comgeom2_mod_h
[5280]	11	USE comdissip_mod_h
[4143]	12	USE infotrac, ONLY: nqtot, tracers
[4056]	13	USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
	14	USE comconst_mod, ONLY: dtvr
	15	USE parallel_lmdz
	16	USE Write_Field_loc
	17	USE Write_Field
	18	USE Bands
	19	USE mod_hallo
	20	USE Vampir
	21	USE times
	22	USE advtrac_mod, ONLY: finmasse
[5271]	23	USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]	24	USE paramet_mod_h
[5271]	25	IMPLICIT NONE
[4056]	26	!
[5271]	27
[5272]	28
[5282]	29	! include "iniprint.h"
[1632]	30
[4056]	31	!---------------------------------------------------------------------------
	32	! Arguments
	33	!---------------------------------------------------------------------------
	34	REAL, INTENT(IN) :: pbarug(ijb_u:ije_u,llm)
	35	REAL, INTENT(IN) :: pbarvg(ijb_v:ije_v,llm)
	36	REAL, INTENT(IN) :: wg(ijb_u:ije_u,llm)
	37	REAL, INTENT(IN) :: p(ijb_u:ije_u,llmp1)
	38	REAL, INTENT(IN) :: massem(ijb_u:ije_u,llm)
	39	REAL, INTENT(INOUT) :: q(ijb_u:ije_u,llm,nqtot)
	40	REAL, INTENT(IN) :: teta(ijb_u:ije_u,llm)
	41	REAL, INTENT(IN) :: pk(ijb_u:ije_u,llm)
	42	!---------------------------------------------------------------------------
	43	! Ajout PPM
	44	!---------------------------------------------------------------------------
	45	REAL :: massebx(ijb_u:ije_u,llm), masseby(ijb_v:ije_v,llm)
	46	!---------------------------------------------------------------------------
	47	! Variables locales
	48	!---------------------------------------------------------------------------
[4064]	49	INTEGER :: ij, l, iq, iadv
[4056]	50	REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu
	51	REAL :: zdp(ijb_u:ije_u), zdpmin, zdpmax
	52	INTEGER, SAVE :: iadvtr=0
	53	!$OMP THREADPRIVATE(iadvtr)
	54	EXTERNAL minmax
[1632]	55
[4056]	56	!---------------------------------------------------------------------------
	57	! Rajouts pour PPM
	58	!---------------------------------------------------------------------------
	59	INTEGER :: indice, n
	60	REAL :: dtbon ! Pas de temps adaptatif pour que CFL<1
	61	REAL :: CFLmaxz, aaa, bbb ! CFL maximum
	62	REAL, DIMENSION(iim,jjb_u:jje_u,llm) :: unatppm, vnatppm, fluxwppm
	63	REAL :: qppm(iim*jjnb_u,llm,nqtot)
	64	REAL :: psppm(iim,jjb_u:jje_u) ! pression au sol
	65	REAL, DIMENSION(llmp1) :: apppm, bpppm
	66	LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE.
	67	INTEGER :: ijb, ije, ijbu, ijbv, ijeu, ijev, j
	68	TYPE(Request),SAVE :: testRequest
[1848]	69	!$OMP THREADPRIVATE(testRequest)
[1632]	70
[4056]	71	! Test sur l'eventuelle creation de valeurs negatives de la masse
	72	ijb = ij_begin; IF(pole_nord) ijb = ij_begin+iip1
	73	ije = ij_end; IF(pole_sud) ije = ij_end-iip1
[1632]	74
[4056]	75	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	76	DO l=1,llm-1
	77	DO ij = ijb+1,ije
	78	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) &
	79	- pbarvg(ij-iip1,l) + pbarvg(ij,l) &
	80	+ wg(ij,l+1) - wg(ij,l)
	81	END DO
	82	! ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
	83	! CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
	84	DO ij = ijb,ije-iip1+1,iip1
	85	zdp(ij)=zdp(ij+iip1-1)
	86	END DO
	87	DO ij = ijb,ije
	88	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
	89	END DO
	90	! CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
	91	! ym ---> eventuellement a revoir
	92	CALL minmax( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
	93	IF(MAX(ABS(zdpmin),ABS(zdpmax)) >0.5) &
	94	WRITE(,)'WARNING DP/P l=',l,' MIN:',zdpmin,' MAX:', zdpmax
	95	END DO
	96	!$OMP END DO NOWAIT
[1632]	97
[4056]	98	!---------------------------------------------------------------------------
	99	! Advection proprement dite (Modification Le Croller (07/2001)
	100	!---------------------------------------------------------------------------
[1632]	101
[4056]	102	!---------------------------------------------------------------------------
	103	! Calcul des moyennes basees sur la masse
	104	!---------------------------------------------------------------------------
	105	!ym CALL massbar_p(massem,massebx,masseby)
[4058]	106	!ym ----> Normalement, inutile pour les schemas classiques
	107	!ym ----> Reverifier lors de la parallelisation des autres schemas
[1632]	108
	109	!
[4056]	110	! CALL Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
	111	! CALL SendRequest(TestRequest)
	112	!!$OMP BARRIER
	113	! CALL WaitRequest(TestRequest)
	114	!$OMP BARRIER
[4058]	115
[4056]	116	! WRITE(,) 'advtrac 157: appel de vlspltgen_loc'
[4058]	117	CALL vlspltgen_loc(q, 2., massem, wg, pbarug, pbarvg, dtvr, p, pk, teta )
[1632]	118
[4056]	119	GOTO 1234
	120	!-------------------------------------------------------------------------
	121	! Appel des sous programmes d'advection
	122	!-------------------------------------------------------------------------
	123	DO iq = 1, nqtot
	124	! CALL clock(t_initial)
	125	IF(tracers(iq)%parent /= 'air') CYCLE
	126	iadv = tracers(iq)%iadv
	127	!-----------------------------------------------------------------------
	128	SELECT CASE(iadv)
	129	!-----------------------------------------------------------------------
	130	CASE(0); CYCLE
	131	!--------------------------------------------------------------------
	132	CASE(10) !--- Schema de Van Leer I MUSCL
	133	!--------------------------------------------------------------------
	134	! WRITE(,) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
	135	!LF CALL vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
[1632]	136
[4056]	137	!--------------------------------------------------------------------
	138	CASE(14) !--- Schema "pseuDO amont" + test sur humidite specifique
	139	!--- pour la vapeur d'eau. F. Codron
	140	!--------------------------------------------------------------------
	141	! WRITE(,) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
[4469]	142	CALL abort_gcm("advtrac","appel a vlspltqs :schema non parallelise",1)
[4056]	143	!LF CALL vlspltqs_p(q(1,1,1),2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta )
[1632]	144
[4056]	145	!--------------------------------------------------------------------
	146	CASE(12) !--- Schema de Frederic Hourdin
	147	!--------------------------------------------------------------------
[4469]	148	CALL abort_gcm("advtrac","appel a vlspltqs :schema non parallelise",1)
[4056]	149	CALL adaptdt(iadv,dtbon,n,pbarug,massem) ! pas de temps adaptatif
	150	IF(n > 1) WRITE(,) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
	151	DO indice=1,n
	152	CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
	153	END DO
[1632]	154
[4056]	155	!--------------------------------------------------------------------
	156	CASE(13) !--- Pas de temps adaptatif
	157	!--------------------------------------------------------------------
[4469]	158	CALL abort_gcm("advtrac","schema non parallelise",1)
[4056]	159	CALL adaptdt(iadv,dtbon,n,pbarug,massem)
	160	IF(n > 1) WRITE(,) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
	161	DO indice=1,n
	162	CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
	163	END DO
[1632]	164
[4056]	165	!--------------------------------------------------------------------
	166	CASE(20) !--- Schema de pente SLOPES
	167	!--------------------------------------------------------------------
[4469]	168	CALL abort_gcm("advtrac","schema SLOPES non parallelise",1)
[4056]	169	CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[1632]	170
[4056]	171	!--------------------------------------------------------------------
	172	CASE(30) !--- Schema de Prather
	173	!--------------------------------------------------------------------
[4469]	174	CALL abort_gcm("advtrac","schema prather non parallelise",1)
[4056]	175	! Pas de temps adaptatif
	176	CALL adaptdt(iadv,dtbon,n,pbarug,massem)
	177	IF(n > 1) WRITE(,) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
	178	CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon)
[1632]	179
[4056]	180	!--------------------------------------------------------------------
	181	CASE(11,16,17,18) !--- Schemas PPM Lin et Rood
	182	!--------------------------------------------------------------------
[4469]	183	CALL abort_gcm("advtrac","schema PPM non parallelise",1)
[4056]	184	! Test sur le flux horizontal
	185	CALL adaptdt(iadv,dtbon,n,pbarug,massem) ! pas de temps adaptatif
	186	IF(n > 1) WRITE(,) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
	187	! Test sur le flux vertical
	188	CFLmaxz=0.
	189	DO l=2,llm
	190	DO ij=iip2,ip1jm
	191	aaa=wg(ij,l)*dtvr/massem(ij,l)
	192	CFLmaxz=max(CFLmaxz,aaa)
	193	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
	194	CFLmaxz=max(CFLmaxz,bbb)
	195	END DO
	196	END DO
	197	IF(CFLmaxz.GE.1) WRITE(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
	198	!----------------------------------------------------------------
	199	! Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
	200	!----------------------------------------------------------------
	201	CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
	202	apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
	203	unatppm,vnatppm,psppm)
[1632]	204
[4056]	205	!----------------------------------------------------------------
	206	DO indice=1,n !--- VL (version PPM) horiz. et PPM vert.
	207	!----------------------------------------------------------------
	208	SELECT CASE(iadv)
	209	!----------------------------------------------------------
	210	CASE(11)
	211	!----------------------------------------------------------
	212	CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
	213	2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
	214	!----------------------------------------------------------
	215	CASE(16) !--- Monotonic PPM
	216	!----------------------------------------------------------
	217	CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
	218	3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
	219	!----------------------------------------------------------
	220	CASE(17) !--- Semi monotonic PPM
	221	!----------------------------------------------------------
	222	CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
	223	4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.)
	224	!----------------------------------------------------------
	225	CASE(18) !--- Positive Definite PPM
	226	!----------------------------------------------------------
	227	CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
	228	5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
	229	END SELECT
	230	!----------------------------------------------------------------
	231	END DO
	232	!----------------------------------------------------------------
	233	! Ss-prg interface PPM3d-LMDZ.4
	234	!----------------------------------------------------------------
	235	CALL interpost(q(1,1,iq),qppm(1,1,iq))
	236	!----------------------------------------------------------------------
	237	END SELECT
	238	!----------------------------------------------------------------------
[1632]	239
[4056]	240	!----------------------------------------------------------------------
[4058]	241	! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
[4056]	242	!----------------------------------------------------------------------
	243	! CALL traceurpole(q(1,1,iq),massem)
[1632]	244
[4056]	245	!--- Calcul du temps cpu pour un schema donne
	246	! CALL clock(t_final)
	247	!ym tps_cpu=t_final-t_initial
	248	!ym cpuadv(iq)=cpuadv(iq)+tps_cpu
[1632]	249
[4056]	250	END DO
[1632]	251
[4056]	252	1234 CONTINUE
	253	!$OMP BARRIER
	254	IF(planet_type=="earth") THEN
[1632]	255	ijb=ij_begin
	256	ije=ij_end
[4056]	257	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	258	DO l = 1, llm
[1632]	259	DO ij = ijb, ije
[4056]	260	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
	261	END DO
	262	END DO
	263	!$OMP END DO
[1632]	264
[4056]	265	CALL qminimum_loc( q, nqtot, finmasse )
[1632]	266
[4056]	267	END IF ! of if (planet_type=="earth")
[1632]	268
[4056]	269	END SUBROUTINE advtrac_loc
[1632]	270

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

Download in other formats:

Original Format