Context Navigation

convadj.F @ 2276

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

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

Download in other formats:

Original Format