Context Navigation

convadj.F @ 3595

Last change on this file since 3595 was 1668, checked in by jvatant, 8 years ago
Another round of cleaning of dust and dummy tracers JVO
File size: 9.1 KB

Rev	Line
[253]	1	subroutine convadj(ngrid,nlay,nq,ptimestep,
	2	& pplay,pplev,ppopsk,
	3	& pu,pv,ph,pq,
	4	& pdufi,pdvfi,pdhfi,pdqfi,
	5	& pduadj,pdvadj,pdhadj,
	6	& pdqadj)
[135]	7
[787]	8	USE tracer_h
[1384]	9	use comcstfi_mod, only: g
[1647]	10	use callkeys_mod, only: tracer
[787]	11
[253]	12	implicit none
[135]	13
[253]	14	!==================================================================
	15	!
	16	! Purpose
	17	! -------
[1647]	18	! Calculates dry convective adjustment.
[253]	19	!
	20	! Authors
	21	! -------
	22	! Original author unknown.
	23	! Modif. 2005 by F. Forget.
	24	!
	25	!==================================================================
[135]	26
[253]	27	! ------------
	28	! Declarations
	29	! ------------
	30
[135]	31
[253]	32	! Arguments
	33	! ---------
[135]	34
	35	INTEGER ngrid,nlay
	36	REAL ptimestep
	37	REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
	38	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
	39	REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
	40	REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)
	41
[253]	42	! Tracers
[135]	43	integer nq
	44	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
	45	real pdqadj(ngrid,nlay,nq)
	46
	47
[253]	48	! Local
	49	! -----
[135]	50
	51	INTEGER ig,i,l,l1,l2,jj
[787]	52	INTEGER jcnt, jadrs(ngrid)
[135]	53
[1308]	54	REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
	55	REAL zu(ngrid,nlay),zv(ngrid,nlay)
	56	REAL zh(ngrid,nlay)
	57	REAL zu2(ngrid,nlay),zv2(ngrid,nlay)
	58	REAL zh2(ngrid,nlay), zhc(ngrid,nlay)
[135]	59	REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
	60
[253]	61	! Tracers
[1647]	62	INTEGER iq
[1308]	63	REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
[1647]	64	REAL zqm(nq)
	65	real A , B
[135]	66	save A, B
[1315]	67	!$OMP THREADPRIVATE(A,B)
[135]	68
	69	real mtot1, mtot2 , mm1, mm2
	70	integer l1ref, l2ref
[1647]	71	LOGICAL vtest(ngrid),down, firstcall
[135]	72	save firstcall
	73	data firstcall/.true./
[1315]	74	!$OMP THREADPRIVATE(firstcall)
[135]	75
[253]	76	! for conservation test
	77	real masse,cadjncons
	78
[135]	79	EXTERNAL SCOPY
[253]	80
	81	! --------------
	82	! Initialisation
	83	! --------------
	84
[135]	85	IF (firstcall) THEN
	86	firstcall=.false.
	87	ENDIF ! of IF (firstcall)
[1647]	88
[135]	89	DO l=1,nlay
	90	DO ig=1,ngrid
	91	zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
	92	zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
	93	zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
	94	ENDDO
	95	ENDDO
	96
	97	if(tracer) then
	98	DO iq =1, nq
	99	DO l=1,nlay
	100	DO ig=1,ngrid
	101	zq(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
	102	ENDDO
	103	ENDDO
	104	ENDDO
	105	end if
	106
	107	CALL scopy(ngrid*nlay,zh,1,zh2,1)
	108	CALL scopy(ngrid*nlay,zu,1,zu2,1)
	109	CALL scopy(ngrid*nlay,zv,1,zv2,1)
	110	CALL scopy(ngridnlaynq,zq,1,zq2,1)
	111
[253]	112	! -----------------------------
	113	! Detection of unstable columns
	114	! -----------------------------
	115	! If ph(above) < ph(below) we set vtest=.true.
	116
[135]	117	DO ig=1,ngrid
[253]	118	vtest(ig)=.false.
[135]	119	ENDDO
[253]	120
[1647]	121	CALL scopy(ngrid*nlay,zh2,1,zhc,1)
[135]	122
[253]	123	! Find out which grid points are convectively unstable
[135]	124	DO l=2,nlay
	125	DO ig=1,ngrid
[253]	126	IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=.true.
[135]	127	ENDDO
	128	ENDDO
[253]	129
	130	! Make a list of them
[135]	131	jcnt=0
	132	DO ig=1,ngrid
	133	IF(vtest(ig)) THEN
	134	jcnt=jcnt+1
	135	jadrs(jcnt)=ig
	136	ENDIF
	137	ENDDO
	138
	139
[253]	140	! ---------------------------------------------------------------
	141	! Adjustment of the "jcnt" unstable profiles indicated by "jadrs"
	142	! ---------------------------------------------------------------
	143
[135]	144	DO jj = 1, jcnt ! loop on every convective grid point
[253]	145
[135]	146	i = jadrs(jj)
	147
[253]	148	! Calculate sigma in this column
[135]	149	DO l=1,nlay+1
	150	sig(l)=pplev(i,l)/pplev(i,1)
	151
	152	ENDDO
	153	DO l=1,nlay
	154	dsig(l)=sig(l)-sig(l+1)
	155	sdsig(l)=ppopsk(i,l)*dsig(l)
	156	ENDDO
	157	l2 = 1
[253]	158
	159	! Test loop upwards on l2
	160
[135]	161	DO
	162	l2 = l2 + 1
	163	IF (l2 .GT. nlay) EXIT
	164	IF (zhc(i, l2) .LT. zhc(i, l2-1)) THEN
	165
[253]	166	! l2 is the highest level of the unstable column
[135]	167
	168	l1 = l2 - 1
	169	l = l1
	170	zsm = sdsig(l2)
	171	zdsm = dsig(l2)
	172	zhm = zh2(i, l2)
[253]	173
	174	! Test loop downwards
	175
[135]	176	DO
	177	zsm = zsm + sdsig(l)
	178	zdsm = zdsm + dsig(l)
	179	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
[1647]	180	zhmc = zhm
[135]	181
[253]	182	! do we have to extend the column downwards?
[135]	183
[253]	184	down = .false.
	185	IF (l1 .ne. 1) then !-- and then
	186	IF (zhmc .lt. zhc(i, l1-1)) then
	187	down = .true.
[135]	188	END IF
	189	END IF
	190
[253]	191	! this could be a problem...
	192
	193	if (down) then
[135]	194
	195	l1 = l1 - 1
	196	l = l1
	197
[253]	198	else
[135]	199
[253]	200	! can we extend the column upwards?
[135]	201
[253]	202	if (l2 .eq. nlay) exit
[135]	203
[253]	204	if (zhc(i, l2+1) .ge. zhmc) exit
	205
[135]	206	l2 = l2 + 1
	207	l = l2
[253]	208
	209	end if
	210
	211	enddo
	212
	213	! New constant profile (average value)
	214
	215
[135]	216	zalpha=0.
	217	zum=0.
	218	zvm=0.
	219	do iq=1,nq
	220	zqm(iq) = 0.
	221	end do
	222	DO l = l1, l2
[1647]	223	zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
[135]	224	zh2(i, l) = zhm
[253]	225	! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	226	zum=zum+dsig(l)*zu2(i,l)
	227	zvm=zvm+dsig(l)*zv2(i,l)
	228	! zum=zum+dsig(l)*zu(i,l)
	229	! zvm=zvm+dsig(l)*zv(i,l)
[135]	230	do iq=1,nq
[253]	231	zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
	232	! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
	233	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	234
	235	! to conserve tracers/ KE, we must calculate zum, zvm and zqm using
	236	! the up-to-date column values. If we do not do this, there are cases
	237	! where convection stops at one level and starts at the next where we
	238	! can break conservation of stuff (particularly tracers) significantly.
	239
[135]	240	end do
	241	ENDDO
	242	zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
	243	zum=zum/(sig(l1)-sig(l2+1))
	244	zvm=zvm/(sig(l1)-sig(l2+1))
	245	do iq=1,nq
	246	zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
	247	end do
	248
	249	IF(zalpha.GT.1.) THEN
	250	zalpha=1.
	251	ELSE
[253]	252	! IF(zalpha.LT.0.) STOP
[135]	253	IF(zalpha.LT.1.e-4) zalpha=1.e-4
	254	ENDIF
	255
	256	DO l=l1,l2
	257	zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
	258	zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
	259	do iq=1,nq
[253]	260	! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))
[135]	261	zq2(i,l,iq)=zqm(iq)
	262	end do
	263	ENDDO
	264
[253]	265
[135]	266	l2 = l2 + 1
[253]	267
	268	END IF ! End of l1 to l2 instability treatment
	269	! We now continue to test from l2 upwards
	270
	271	ENDDO ! End of upwards loop on l2
	272
	273
	274	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	275	! check conservation
	276	cadjncons=0.0
	277
	278	if(cadjncons.lt.-1.e-6)then
	279	print*,'convadj has just crashed...'
	280	print*,'i = ',i
	281	print*,'l1 = ',l1
	282	print*,'l2 = ',l2
	283	print*,'cadjncons = ',cadjncons
	284	do l = 1, nlay
	285	print*,'dsig = ',dsig(l)
	286	end do
	287	do l = 1, nlay
	288	print*,'dsig = ',dsig(l)
	289	end do
	290	do l = 1, nlay
	291	print*,'sig = ',sig(l)
	292	end do
	293	do l = 1, nlay
	294	print*,'pplay(ig,:) = ',pplay(i,l)
	295	end do
	296	do l = 1, nlay+1
	297	print*,'pplev(ig,:) = ',pplev(i,l)
	298	end do
	299	do l = 1, nlay
	300	print*,'ph(ig,:) = ',ph(i,l)
	301	end do
	302	do l = 1, nlay
	303	print*,'ph(ig,:) = ',ph(i,l)
	304	end do
	305	do l = 1, nlay
	306	print*,'ph(ig,:) = ',ph(i,l)
	307	end do
	308	do l = 1, nlay
	309	print*,'zh(ig,:) = ',zh(i,l)
	310	end do
	311	do l = 1, nlay
	312	print*,'zh2(ig,:) = ',zh2(i,l)
	313	end do
[1668]	314	print*,'jadrs=',jadrs
[253]	315
	316	call abort
	317	endif
	318	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	319
	320
[135]	321	ENDDO
[253]	322
[135]	323	DO l=1,nlay
	324	DO ig=1,ngrid
	325	pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
	326	pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
	327	pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
	328	ENDDO
	329	ENDDO
	330
	331	if(tracer) then
	332	do iq=1, nq
	333	do l=1,nlay
[787]	334	DO ig=1,ngrid
[135]	335	pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep
	336	end do
	337	end do
	338	end do
	339	end if
	340
[253]	341
	342	! output
[135]	343	! if (ngrid.eq.1) then
	344	! ig=1
	345	! iq =1
	346	! write(,)'**** l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)'
	347	! do l=nlay,1,-1
	348	! write(,) l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)
	349	! end do
	350	! end if
	351
	352
[253]	353	return
	354	end

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

Download in other formats:

Original Format