Context Navigation

convadj.F @ 735

Last change on this file since 735 was 253, checked in by emillour, 14 years ago

Generic GCM

Massive update to version 0.7

EM+RW

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

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

Download in other formats:

Original Format