Context Navigation

convadj.F @ 3276

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

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

Download in other formats:

Original Format