Context Navigation

integrd_p.F @ 2902

Last change on this file since 2902 was 2126, checked in by emillour, 6 years ago

Common dynamics:
Some work to enforce total tracer mass conservation in the dynamics.
Still to be further studied and validated.
For now these changes are triggered by setting a "force_conserv_tracer"
flag to ".true." in run.def (default is ".false." to not change anything
with respect to previous versions).
When force_conserv_tracer=.true. then:

Rescale tracer mass in caladvtrac after tracer advection computations
Recompute q ratios once atmospheric mass has been updated in integrd

These steps technically ensure total tracer mass conservation but it
might be the tracer advection scheme and/or time-stepping updating
sequence of fields that should be rethought or fixed.
EM

File size: 10.4 KB

Rev	Line
[1]	1	!
[776]	2	! $Id: integrd_p.F 1616 2012-02-17 11:59:00Z emillour $
[1]	3	!
	4	SUBROUTINE integrd_p
	5	$ ( nq,vcovm1,ucovm1,tetam1,psm1,massem1,
[776]	6	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis) !,finvmaold)
[1510]	7	USE parallel_lmdz, ONLY: ij_begin, ij_end, pole_nord, pole_sud,
	8	& omp_chunk
[2126]	9	USE control_mod, only : planet_type,force_conserv_tracer
[1422]	10	USE comvert_mod, ONLY: ap,bp
	11	USE comconst_mod, ONLY: pi
	12	USE logic_mod, ONLY: leapf
	13	USE temps_mod, ONLY: dt
[1]	14	IMPLICIT NONE
	15
	16
	17	c=======================================================================
	18	c
	19	c Auteur: P. Le Van
	20	c -------
	21	c
	22	c objet:
	23	c ------
	24	c
	25	c Incrementation des tendances dynamiques
	26	c
	27	c=======================================================================
	28	c-----------------------------------------------------------------------
	29	c Declarations:
	30	c -------------
	31
	32	#include "dimensions.h"
	33	#include "paramet.h"
	34	#include "comgeom.h"
[776]	35	#include "iniprint.h"
[1]	36
	37	c Arguments:
	38	c ----------
	39
[776]	40	integer,intent(in) :: nq ! number of tracers to handle in this routine
	41	real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind
	42	real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind
	43	real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature
	44	real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers
	45	real,intent(inout) :: ps0(ip1jmp1) ! surface pressure
	46	real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass
	47	real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
	48	! values at previous time step
	49	real,intent(inout) :: vcovm1(ip1jm,llm)
	50	real,intent(inout) :: ucovm1(ip1jmp1,llm)
	51	real,intent(inout) :: tetam1(ip1jmp1,llm)
	52	real,intent(inout) :: psm1(ip1jmp1)
	53	real,intent(inout) :: massem1(ip1jmp1,llm)
	54	! the tendencies to add
	55	real,intent(in) :: dv(ip1jm,llm)
	56	real,intent(in) :: du(ip1jmp1,llm)
	57	real,intent(in) :: dteta(ip1jmp1,llm)
	58	real,intent(in) :: dp(ip1jmp1)
	59	real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused
	60	! real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused
[1]	61
	62	c Local:
	63	c ------
	64
	65	REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
[776]	66	REAL massescr( ip1jmp1,llm )
[2126]	67	REAL :: massratio(ip1jmp1,llm)
[776]	68	! REAL finvmasse(ip1jmp1,llm)
[1]	69	REAL,SAVE :: p(ip1jmp1,llmp1)
	70	REAL tpn,tps,tppn(iim),tpps(iim)
	71	REAL qpn,qps,qppn(iim),qpps(iim)
	72	REAL,SAVE :: deltap( ip1jmp1,llm )
	73
[776]	74	INTEGER l,ij,iq,i,j
[1]	75
	76	REAL SSUM
	77	EXTERNAL SSUM
	78	INTEGER ijb,ije,jjb,jje
[1959]	79	REAL,SAVE :: ps(ip1jmp1)=0
[1]	80	LOGICAL :: checksum
	81	INTEGER :: stop_it
	82	c-----------------------------------------------------------------------
	83	c$OMP BARRIER
	84	if (pole_nord) THEN
	85	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	86	DO l = 1,llm
	87	DO ij = 1,iip1
	88	ucov( ij , l) = 0.
	89	uscr( ij ) = 0.
	90	ENDDO
	91	ENDDO
	92	c$OMP END DO NOWAIT
	93	ENDIF
	94
	95	if (pole_sud) THEN
	96	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	97	DO l = 1,llm
	98	DO ij = 1,iip1
	99	ucov( ij +ip1jm, l) = 0.
	100	uscr( ij +ip1jm ) = 0.
	101	ENDDO
	102	ENDDO
	103	c$OMP END DO NOWAIT
	104	ENDIF
	105
	106	c ............ integration de ps ..............
	107
	108	c CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
	109
	110	ijb=ij_begin
	111	ije=ij_end
	112	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	113	DO l = 1,llm
	114	massescr(ijb:ije,l)=masse(ijb:ije,l)
	115	ENDDO
	116	c$OMP END DO NOWAIT
	117
	118	c$OMP DO SCHEDULE(STATIC)
	119	DO 2 ij = ijb,ije
	120	pscr (ij) = ps0(ij)
	121	ps (ij) = psm1(ij) + dt * dp(ij)
	122	2 CONTINUE
	123	c$OMP END DO
	124	c$OMP BARRIER
	125	c --> ici synchro OPENMP pour ps
	126
	127	checksum=.TRUE.
	128	stop_it=0
	129
	130	c$OMP DO SCHEDULE(STATIC)
	131	DO ij = ijb,ije
	132	IF( ps(ij).LT.0. ) THEN
	133	IF (checksum) stop_it=ij
	134	checksum=.FALSE.
	135	ENDIF
	136	ENDDO
	137	c$OMP END DO NOWAIT
	138
	139	IF( .NOT. checksum ) THEN
[776]	140	write(lunout,*) "integrd: negative surface pressure ",
	141	& ps(stop_it)
	142	write(lunout,*) " at node ij =", stop_it
	143	! since ij=j+(i-1)*jjp1 , we have
	144	j=modulo(stop_it,jjp1)
	145	i=1+(stop_it-j)/jjp1
	146	write(lunout,) " lon = ",rlonv(i)180./pi, " deg",
	147	& " lat = ",rlatu(j)*180./pi, " deg"
[907]	148	write(lunout,*) " psm1(ij)=",psm1(stop_it)," dt=",dt,
	149	& " dp(ij)=",dp(stop_it)
[1391]	150	call abort_gcm("integrd_p", "negative surface pressure", 1)
[1]	151	ENDIF
	152
	153	c
	154	C$OMP MASTER
	155	if (pole_nord) THEN
	156
	157	DO ij = 1, iim
	158	tppn(ij) = aire( ij ) * ps( ij )
	159	ENDDO
	160	tpn = SSUM(iim,tppn,1)/apoln
	161	DO ij = 1, iip1
	162	ps( ij ) = tpn
	163	ENDDO
	164
	165	ENDIF
	166
	167	if (pole_sud) THEN
	168
	169	DO ij = 1, iim
	170	tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
	171	ENDDO
	172	tps = SSUM(iim,tpps,1)/apols
	173	DO ij = 1, iip1
	174	ps(ij+ip1jm) = tps
	175	ENDDO
	176
	177	ENDIF
	178	c$OMP END MASTER
	179	c$OMP BARRIER
	180	c
	181	c ... Calcul de la nouvelle masse d'air au dernier temps integre t+1 ...
	182	c
	183
	184	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
	185	c$OMP BARRIER
	186	CALL massdair_p ( p , masse )
	187
[776]	188	! Ehouarn : we don't use/need finvmaold and finvmasse,
	189	! so might as well not compute them
	190	!c CALL SCOPY( ijp1llm , masse, 1, finvmasse, 1 )
	191	! ijb=ij_begin
	192	! ije=ij_end
	193	!
	194	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	195	! DO l = 1,llm
	196	! finvmasse(ijb:ije,l)=masse(ijb:ije,l)
	197	! ENDDO
	198	!c$OMP END DO NOWAIT
	199	!
	200	! jjb=jj_begin
	201	! jje=jj_end
	202	! CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1 )
[1]	203	c
	204
	205	c ............ integration de ucov, vcov, h ..............
	206
	207	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	208	DO 10 l = 1,llm
	209
	210	ijb=ij_begin
	211	ije=ij_end
	212	if (pole_nord) ijb=ij_begin+iip1
	213	if (pole_sud) ije=ij_end-iip1
	214
	215	DO 4 ij = ijb,ije
	216	uscr( ij ) = ucov( ij,l )
	217	ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
	218	4 CONTINUE
	219
	220	ijb=ij_begin
	221	ije=ij_end
	222	if (pole_sud) ije=ij_end-iip1
	223
	224	DO 5 ij = ijb,ije
	225	vscr( ij ) = vcov( ij,l )
	226	vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
	227	5 CONTINUE
	228
	229	ijb=ij_begin
	230	ije=ij_end
	231
	232	DO 6 ij = ijb,ije
	233	hscr( ij ) = teta(ij,l)
	234	teta ( ij,l ) = tetam1(ij,l) * massem1(ij,l) / masse(ij,l)
	235	$ + dt * dteta(ij,l) / masse(ij,l)
	236	6 CONTINUE
	237
	238	c .... Calcul de la valeur moyenne, unique aux poles pour teta ......
	239	c
	240	c
	241	IF (pole_nord) THEN
	242
	243	DO ij = 1, iim
	244	tppn(ij) = aire( ij ) * teta( ij ,l)
	245	ENDDO
	246	tpn = SSUM(iim,tppn,1)/apoln
	247
	248	DO ij = 1, iip1
	249	teta( ij ,l) = tpn
	250	ENDDO
	251
	252	ENDIF
	253
	254	IF (pole_sud) THEN
	255
	256	DO ij = 1, iim
	257	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	258	ENDDO
	259	tps = SSUM(iim,tpps,1)/apols
	260
	261	DO ij = 1, iip1
	262	teta(ij+ip1jm,l) = tps
	263	ENDDO
	264
	265	ENDIF
	266	c
	267
	268	IF(leapf) THEN
	269	c CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
	270	c CALL SCOPY ( ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
	271	c CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
	272	ijb=ij_begin
	273	ije=ij_end
	274	ucovm1(ijb:ije,l)=uscr(ijb:ije)
	275	tetam1(ijb:ije,l)=hscr(ijb:ije)
	276	if (pole_sud) ije=ij_end-iip1
	277	vcovm1(ijb:ije,l)=vscr(ijb:ije)
	278
	279	END IF
	280
	281	10 CONTINUE
	282	c$OMP END DO NOWAIT
	283
	284	c
	285	c ....... integration de q ......
	286	c
	287	ijb=ij_begin
	288	ije=ij_end
	289
	290	if (planet_type.eq."earth") then
	291	! Earth-specific treatment of first 2 tracers (water)
	292	c$OMP BARRIER
	293	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	294	DO l = 1, llm
	295	DO ij = ijb, ije
	296	deltap(ij,l) = p(ij,l) - p(ij,l+1)
	297	ENDDO
	298	ENDDO
	299	c$OMP END DO NOWAIT
	300	c$OMP BARRIER
	301
	302	CALL qminimum_p( q, nq, deltap )
	303	c
	304	c ..... Calcul de la valeur moyenne, unique aux poles pour q .....
	305	c
	306	c$OMP BARRIER
	307	IF (pole_nord) THEN
	308
	309	DO iq = 1, nq
	310
	311	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	312	DO l = 1, llm
	313
	314	DO ij = 1, iim
	315	qppn(ij) = aire( ij ) * q( ij ,l,iq)
	316	ENDDO
	317	qpn = SSUM(iim,qppn,1)/apoln
	318
	319	DO ij = 1, iip1
	320	q( ij ,l,iq) = qpn
	321	ENDDO
	322
	323	ENDDO
	324	c$OMP END DO NOWAIT
	325
	326	ENDDO
	327
	328	ENDIF
	329
	330	IF (pole_sud) THEN
	331
	332	DO iq = 1, nq
	333
	334	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	335	DO l = 1, llm
	336
	337	DO ij = 1, iim
	338	qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
	339	ENDDO
	340	qps = SSUM(iim,qpps,1)/apols
	341
	342	DO ij = 1, iip1
	343	q(ij+ip1jm,l,iq) = qps
	344	ENDDO
	345
	346	ENDDO
	347	c$OMP END DO NOWAIT
	348
	349	ENDDO
	350
	351	ENDIF
	352
[776]	353	! Ehouarn: forget about finvmaold
	354	!c CALL SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
	355	!
	356	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	357	! DO l = 1, llm
	358	! finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)
	359	! ENDDO
	360	!c$OMP END DO NOWAIT
[1]	361
[7]	362	endif ! of if (planet_type.eq."earth")
	363
[2126]	364
	365	if (force_conserv_tracer) then
	366	! Ehouarn: try to keep total amont of tracers fixed
	367	! by acounting for mass change in each cell
	368	massratio(ijb:ije,1:llm)=massescr(ijb:ije,1:llm)
	369	& /masse(ijb:ije,1:llm)
	370	do iq=1,nq
	371	q(ijb:ije,1:llm,iq)=q(ijb:ije,1:llm,iq)
	372	& *massratio(ijb:ije,1:llm)
	373	enddo
	374	endif ! of if (force_conserv_tracer)
[1]	375	c
	376	c
	377	c ..... FIN de l'integration de q .......
	378
	379	15 continue
	380
	381	c$OMP DO SCHEDULE(STATIC)
	382	DO ij=ijb,ije
	383	ps0(ij)=ps(ij)
	384	ENDDO
	385	c$OMP END DO NOWAIT
	386
	387	c .................................................................
	388
	389
	390	IF( leapf ) THEN
	391	c CALL SCOPY ( ip1jmp1 , pscr , 1, psm1 , 1 )
	392	c CALL SCOPY ( ip1jmp1*llm, massescr, 1, massem1, 1 )
	393	c$OMP DO SCHEDULE(STATIC)
	394	DO ij=ijb,ije
	395	psm1(ij)=pscr(ij)
	396	ENDDO
	397	c$OMP END DO NOWAIT
	398
	399	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	400	DO l = 1, llm
	401	massem1(ijb:ije,l)=massescr(ijb:ije,l)
	402	ENDDO
	403	c$OMP END DO NOWAIT
	404	END IF
	405	c$OMP BARRIER
	406	RETURN
	407	END

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

Download in other formats:

Original Format