Context Navigation

convadj.F @ 1364

Last change on this file since 1364 was 1315, checked in by milmd, 12 years ago
LMDZ.GENERIC. OpenMP directives added in generic physic. When running in pure OpenMP or hybrid OpenMP/MPI, may have some bugs with condense_cloud and wstats routines.
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,
	5	& pduadj,pdvadj,pdhadj,
	6	& pdqadj)
[135]	7
[787]	8	USE tracer_h
	9
[253]	10	implicit none
[135]	11
[253]	12	!==================================================================
	13	!
	14	! Purpose
	15	! -------
	16	! Calculates dry convective adjustment. If one tracer is CO2,
	17	! we take into account the molecular mass variation (e.g. when
	18	! CO2 condenses) to trigger convection (F. Forget 01/2005)
	19	!
	20	! Authors
	21	! -------
	22	! Original author unknown.
	23	! Modif. 2005 by F. Forget.
	24	!
	25	!==================================================================
[135]	26
[253]	27	! ------------
	28	! Declarations
	29	! ------------
	30
[1308]	31	!#include "dimensions.h"
	32	!#include "dimphys.h"
[135]	33	#include "comcstfi.h"
	34	#include "callkeys.h"
	35
	36
[253]	37	! Arguments
	38	! ---------
[135]	39
	40	INTEGER ngrid,nlay
	41	REAL ptimestep
	42	REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
	43	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
	44	REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
	45	REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)
	46
[253]	47	! Tracers
[135]	48	integer nq
	49	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
	50	real pdqadj(ngrid,nlay,nq)
	51
	52
[253]	53	! Local
	54	! -----
[135]	55
	56	INTEGER ig,i,l,l1,l2,jj
[787]	57	INTEGER jcnt, jadrs(ngrid)
[135]	58
[1308]	59	REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
	60	REAL zu(ngrid,nlay),zv(ngrid,nlay)
	61	REAL zh(ngrid,nlay)
	62	REAL zu2(ngrid,nlay),zv2(ngrid,nlay)
	63	REAL zh2(ngrid,nlay), zhc(ngrid,nlay)
[135]	64	REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
	65
[253]	66	! Tracers
[135]	67	INTEGER iq,ico2
	68	save ico2
[1315]	69	!$OMP THREADPRIVATE(ico2)
[1308]	70	REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
[787]	71	REAL zqm(nq),zqco2m
[135]	72	real m_co2, m_noco2, A , B
	73	save A, B
[1315]	74	!$OMP THREADPRIVATE(A,B)
[135]	75
	76	real mtot1, mtot2 , mm1, mm2
	77	integer l1ref, l2ref
[787]	78	LOGICAL vtest(ngrid),down,firstcall
[135]	79	save firstcall
	80	data firstcall/.true./
[1315]	81	!$OMP THREADPRIVATE(firstcall)
[135]	82
[253]	83	! for conservation test
	84	real masse,cadjncons
	85
[135]	86	EXTERNAL SCOPY
[253]	87
	88	! --------------
	89	! Initialisation
	90	! --------------
	91
[135]	92	IF (firstcall) THEN
	93	ico2=0
	94	if (tracer) then
[253]	95	! Prepare Special treatment if one of the tracers is CO2 gas
[787]	96	do iq=1,nq
[135]	97	if (noms(iq).eq."co2") then
	98	print*,'dont go there'
	99	stop
	100	ico2=iq
	101	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
	102	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
[253]	103	! Compute A and B coefficient use to compute
	104	! mean molecular mass Mair defined by
	105	! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
	106	! 1/Mair = A*q(ico2) + B
[135]	107	A =(1/m_co2 - 1/m_noco2)
	108	B=1/m_noco2
	109	end if
	110	enddo
	111	endif
	112	firstcall=.false.
	113	ENDIF ! of IF (firstcall)
	114
	115	DO l=1,nlay
	116	DO ig=1,ngrid
	117	zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
	118	zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
	119	zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
	120	ENDDO
	121	ENDDO
	122
	123	if(tracer) then
	124	DO iq =1, nq
	125	DO l=1,nlay
	126	DO ig=1,ngrid
	127	zq(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
	128	ENDDO
	129	ENDDO
	130	ENDDO
	131	end if
	132
	133	CALL scopy(ngrid*nlay,zh,1,zh2,1)
	134	CALL scopy(ngrid*nlay,zu,1,zu2,1)
	135	CALL scopy(ngrid*nlay,zv,1,zv2,1)
	136	CALL scopy(ngridnlaynq,zq,1,zq2,1)
	137
[253]	138	! -----------------------------
	139	! Detection of unstable columns
	140	! -----------------------------
	141	! If ph(above) < ph(below) we set vtest=.true.
	142
[135]	143	DO ig=1,ngrid
[253]	144	vtest(ig)=.false.
[135]	145	ENDDO
[253]	146
[135]	147	if (ico2.ne.0) then
[253]	148	! Special case if one of the tracers is CO2 gas
[135]	149	DO l=1,nlay
	150	DO ig=1,ngrid
	151	zhc(ig,l) = zh2(ig,l)(Azq2(ig,l,ico2)+B)
	152	ENDDO
	153	ENDDO
	154	else
	155	CALL scopy(ngrid*nlay,zh2,1,zhc,1)
	156	end if
	157
[253]	158	! Find out which grid points are convectively unstable
[135]	159	DO l=2,nlay
	160	DO ig=1,ngrid
[253]	161	IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=.true.
[135]	162	ENDDO
	163	ENDDO
[253]	164
	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
	199	IF (zhc(i, l2) .LT. zhc(i, l2-1)) THEN
	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
	228	IF (zhmc .lt. zhc(i, l1-1)) then
	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