Context Navigation

convadj.F @ 3506

Last change on this file since 3506 was 3275, checked in by afalco, 8 months ago
Pluto PCM: Changed _vap to _gas; Included surfprop.F90; callcorrk includes methane AF
File size: 9.2 KB

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

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

Download in other formats:

Original Format