Context Navigation

convadj.F @ 2176

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

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

Download in other formats:

Original Format