Context Navigation

convadj.F @ 3469

Last change on this file since 3469 was 3322, checked in by emillour, 8 months ago
Generic PCM: Cleanup and cosmetics on convective adjustment routine. Make it a module in the process. EM
File size: 10.7 KB

Rev	Line
[3322]	1	module convadj_mod
	2
	3	implicit none
	4
	5	contains
	6
[253]	7	subroutine convadj(ngrid,nlay,nq,ptimestep,
	8	& pplay,pplev,ppopsk,
	9	& pu,pv,ph,pq,
	10	& pdufi,pdvfi,pdhfi,pdqfi,
[2232]	11	& pduadj,pdvadj,pdhadj,pdqadj)
[135]	12
[3322]	13	use tracer_h, only: igcm_h2o_vap
[1384]	14	use comcstfi_mod, only: g
[3322]	15	use generic_cloud_common_h, only: epsi_generic
[3236]	16	use generic_tracer_index_mod, only: generic_tracer_index
	17	use callkeys_mod, only: tracer,water,generic_condensation,
	18	& virtual_correction
[787]	19
[253]	20	implicit none
[135]	21
[253]	22	!==================================================================
	23	!
	24	! Purpose
	25	! -------
[3322]	26	! Compute dry convective adjustment.
	27	! See old reference paper: Hourdin et al. JAS 1993
	28	! "Meteorological Variability and the Annual Surface
	29	! Pressure Cycle on Mars"
	30	! https://doi.org/10.1175/1520-0469(1993)050%3C3625:MVATAS%3E2.0.CO;2
[253]	31	!==================================================================
[135]	32
[253]	33	! Arguments
	34	! ---------
[135]	35
[3322]	36	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
	37	INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers
	38	INTEGER,INTENT(IN) :: nq ! number of tracers
	39	REAL,INTENT(IN) :: ptimestep ! physics time step (s)
	40	REAL,INTENT(IN) :: pplay(ngrid,nlay) ! mi-layer pressure (Pa)
	41	REAL,INTENT(IN) :: pplev(ngrid,nlay+1) ! inter-layer pressure (Pa)
	42	REAL,INTENT(IN) :: ppopsk(ngrid,nlay) ! Exner (wrt surface pressure)
	43	! fields from dynamics
	44	REAL,INTENT(IN) :: ph(ngrid,nlay) ! potential temperature (K)
	45	REAL,INTENT(IN) :: pu(ngrid,nlay) ! zonal wind (m/s)
	46	REAL,INTENT(IN) :: pv(ngrid,nlay) ! meridioanl wind (m/s)
	47	REAL,INTENT(IN) :: pq(ngrid,nlay,nq) ! tracers (../kg_air)
	48	! tendencies due to previous physical processes
	49	REAL,INTENT(IN) :: pdufi(ngrid,nlay) ! on zonal wind (m/s/s)
	50	REAL,INTENT(IN) :: pdvfi(ngrid,nlay) ! on meridional wind (m/s/s)
	51	REAL,INTENT(IN) :: pdhfi(ngrid,nlay)! on potential temperature (/K/s)
	52	REAL,INTENT(IN) :: pdqfi(ngrid,nlay,nq) ! on tracers (../kg_air/s)
	53	! tendencies due to convetive adjustement
	54	REAL,INTENT(OUT) :: pduadj(ngrid,nlay) ! on zonal wind (m/s/s)
	55	REAL,INTENT(OUT) :: pdvadj(ngrid,nlay) ! on meridinal wind (m/s/s)
	56	REAL,INTENT(OUT) :: pdhadj(ngrid,nlay) ! on potential temperature (/K/s)
	57	REAL,INTENT(OUT) :: pdqadj(ngrid,nlay,nq) ! on traceurs (../kg_air/s)
[135]	58
	59
[253]	60	! Local
	61	! -----
[135]	62
	63	INTEGER ig,i,l,l1,l2,jj
[787]	64	INTEGER jcnt, jadrs(ngrid)
[135]	65
[1308]	66	REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
	67	REAL zu(ngrid,nlay),zv(ngrid,nlay)
[3236]	68	REAL zh(ngrid,nlay), zvh(ngrid,nlay)
[1308]	69	REAL zu2(ngrid,nlay),zv2(ngrid,nlay)
[3236]	70	REAL zh2(ngrid,nlay),zvh2(ngrid,nlay),zhc(ngrid,nlay)
[135]	71	REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
	72
[253]	73	! Tracers
[2482]	74	INTEGER iq
[1308]	75	REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
[3322]	76	REAL zqm(nq)
[3236]	77
	78	integer igcm_generic_vap, igcm_generic_ice
	79	logical call_ice_vap_generic
[135]	80
[2482]	81	LOGICAL vtest(ngrid),down
[135]	82
[253]	83	! for conservation test
	84	real masse,cadjncons
	85
	86	! --------------
	87	! Initialisation
	88	! --------------
[3236]	89	zh(:,:)=ph(:,:)+pdhfi(:,:)*ptimestep
	90	zu(:,:)=pu(:,:)+pdufi(:,:)*ptimestep
	91	zv(:,:)=pv(:,:)+pdvfi(:,:)*ptimestep
[253]	92
[135]	93	if(tracer) then
[3236]	94	zq(:,:,:)=pq(:,:,:)+pdqfi(:,:,:)*ptimestep
[135]	95	end if
	96
[3322]	97	zh2(:,:)=zh(:,:)
	98	zu2(:,:)=zu(:,:)
	99	zv2(:,:)=zv(:,:)
	100	zq2(:,:,:)=zq(:,:,:)
[135]	101
[253]	102	! -----------------------------
	103	! Detection of unstable columns
	104	! -----------------------------
	105	! If ph(above) < ph(below) we set vtest=.true.
	106
[135]	107	DO ig=1,ngrid
[253]	108	vtest(ig)=.false.
[135]	109	ENDDO
[253]	110
[3236]	111	if((generic_condensation) .and. (virtual_correction)) THEN
	112	DO iq=1,nq
	113	call generic_tracer_index(nq,iq,igcm_generic_vap,
	114	& igcm_generic_ice,call_ice_vap_generic)
	115	if(call_ice_vap_generic) then
	116	zvh(:,:)=zh(:,:)*
	117	& (1.+zq(:,:,igcm_generic_vap)/epsi_generic)/
	118	& (1.+zq(:,:,igcm_generic_vap))
	119	endif
	120	ENDDO
[3322]	121	zvh2(:,:)=zvh(:,:)
	122	zhc(:,:)=zvh2(:,:)
[3236]	123	else
[3322]	124	zhc(:,:)=zh2(:,:)
[3236]	125	endif
[135]	126
[253]	127	! Find out which grid points are convectively unstable
[135]	128	DO l=2,nlay
	129	DO ig=1,ngrid
[2232]	130	IF (zhc(ig,l).LT.zhc(ig,l-1)) THEN
[2107]	131	vtest(ig)=.true.
	132	ENDIF
[135]	133	ENDDO
	134	ENDDO
[2107]	135
[253]	136	! Make a list of them
[135]	137	jcnt=0
	138	DO ig=1,ngrid
	139	IF(vtest(ig)) THEN
	140	jcnt=jcnt+1
	141	jadrs(jcnt)=ig
	142	ENDIF
	143	ENDDO
	144
	145
[253]	146	! ---------------------------------------------------------------
	147	! Adjustment of the "jcnt" unstable profiles indicated by "jadrs"
	148	! ---------------------------------------------------------------
	149
[135]	150	DO jj = 1, jcnt ! loop on every convective grid point
[253]	151
[135]	152	i = jadrs(jj)
	153
[253]	154	! Calculate sigma in this column
[135]	155	DO l=1,nlay+1
	156	sig(l)=pplev(i,l)/pplev(i,1)
	157
	158	ENDDO
	159	DO l=1,nlay
	160	dsig(l)=sig(l)-sig(l+1)
	161	sdsig(l)=ppopsk(i,l)*dsig(l)
	162	ENDDO
	163	l2 = 1
[253]	164
	165	! Test loop upwards on l2
	166
[135]	167	DO
	168	l2 = l2 + 1
	169	IF (l2 .GT. nlay) EXIT
[2232]	170	IF (zhc(i, l2).LT.zhc(i, l2-1)) THEN
[135]	171
[253]	172	! l2 is the highest level of the unstable column
[135]	173
	174	l1 = l2 - 1
	175	l = l1
	176	zsm = sdsig(l2)
	177	zdsm = dsig(l2)
	178	zhm = zh2(i, l2)
[253]	179
	180	! Test loop downwards
	181
[135]	182	DO
	183	zsm = zsm + sdsig(l)
	184	zdsm = zdsm + dsig(l)
	185	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
[3236]	186
	187	if (generic_condensation .and. virtual_correction) then
	188	zhmc = zhm*
	189	& (1.+zq(i,l,igcm_generic_vap)/epsi_generic)/
	190	& (1.+zq(i,l,igcm_generic_vap))
	191	else
	192	zhmc = zhm
	193	endif
[135]	194
[253]	195	! do we have to extend the column downwards?
[135]	196
[253]	197	down = .false.
	198	IF (l1 .ne. 1) then !-- and then
[2232]	199	IF (zhmc.LT.zhc(i, l1-1)) then
[253]	200	down = .true.
[135]	201	END IF
	202	END IF
	203
[253]	204	! this could be a problem...
	205
	206	if (down) then
[135]	207
	208	l1 = l1 - 1
	209	l = l1
	210
[253]	211	else
[135]	212
[253]	213	! can we extend the column upwards?
[135]	214
[253]	215	if (l2 .eq. nlay) exit
[135]	216
[253]	217	if (zhc(i, l2+1) .ge. zhmc) exit
	218
[135]	219	l2 = l2 + 1
	220	l = l2
[253]	221
	222	end if
	223
	224	enddo
	225
	226	! New constant profile (average value)
	227
	228
[135]	229	zalpha=0.
	230	zum=0.
	231	zvm=0.
	232	do iq=1,nq
	233	zqm(iq) = 0.
	234	end do
	235	DO l = l1, l2
[2482]	236	zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
[135]	237	zh2(i, l) = zhm
[253]	238	! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	239	zum=zum+dsig(l)*zu2(i,l)
	240	zvm=zvm+dsig(l)*zv2(i,l)
	241	! zum=zum+dsig(l)*zu(i,l)
	242	! zvm=zvm+dsig(l)*zv(i,l)
[135]	243	do iq=1,nq
[253]	244	zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
	245	! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
	246	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	247
	248	! to conserve tracers/ KE, we must calculate zum, zvm and zqm using
	249	! the up-to-date column values. If we do not do this, there are cases
	250	! where convection stops at one level and starts at the next where we
	251	! can break conservation of stuff (particularly tracers) significantly.
	252
[135]	253	end do
	254	ENDDO
	255	zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
	256	zum=zum/(sig(l1)-sig(l2+1))
	257	zvm=zvm/(sig(l1)-sig(l2+1))
	258	do iq=1,nq
	259	zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
	260	end do
	261
	262	IF(zalpha.GT.1.) THEN
	263	zalpha=1.
	264	ELSE
[253]	265	! IF(zalpha.LT.0.) STOP
[135]	266	IF(zalpha.LT.1.e-4) zalpha=1.e-4
	267	ENDIF
	268
	269	DO l=l1,l2
	270	zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
	271	zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
	272	do iq=1,nq
[253]	273	! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))
[135]	274	zq2(i,l,iq)=zqm(iq)
	275	end do
	276	ENDDO
	277
[253]	278
[135]	279	l2 = l2 + 1
[253]	280
	281	END IF ! End of l1 to l2 instability treatment
	282	! We now continue to test from l2 upwards
	283
	284	ENDDO ! End of upwards loop on l2
	285
	286
	287	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	288	! check conservation
	289	cadjncons=0.0
	290	if(water)then
[3322]	291	do l = 1, nlay
[253]	292	masse = (pplev(i,l) - pplev(i,l+1))/g
	293	iq = igcm_h2o_vap
	294	cadjncons = cadjncons +
	295	& masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep
[3322]	296	end do
[253]	297	endif
	298
	299	if(cadjncons.lt.-1.e-6)then
	300	print*,'convadj has just crashed...'
	301	print*,'i = ',i
	302	print*,'l1 = ',l1
	303	print*,'l2 = ',l2
	304	print*,'cadjncons = ',cadjncons
	305	do l = 1, nlay
	306	print*,'dsig = ',dsig(l)
	307	end do
	308	do l = 1, nlay
	309	print*,'dsig = ',dsig(l)
	310	end do
	311	do l = 1, nlay
	312	print*,'sig = ',sig(l)
	313	end do
	314	do l = 1, nlay
	315	print*,'pplay(ig,:) = ',pplay(i,l)
	316	end do
	317	do l = 1, nlay+1
	318	print*,'pplev(ig,:) = ',pplev(i,l)
	319	end do
	320	do l = 1, nlay
	321	print*,'ph(ig,:) = ',ph(i,l)
	322	end do
	323	do l = 1, nlay
	324	print*,'ph(ig,:) = ',ph(i,l)
	325	end do
	326	do l = 1, nlay
	327	print*,'ph(ig,:) = ',ph(i,l)
	328	end do
	329	do l = 1, nlay
	330	print*,'zh(ig,:) = ',zh(i,l)
	331	end do
	332	do l = 1, nlay
	333	print*,'zh2(ig,:) = ',zh2(i,l)
	334	end do
	335	do l = 1, nlay
	336	print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_vap)
	337	end do
	338	do l = 1, nlay
	339	print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_vap)
	340	end do
	341	print*,'zqm(vap) = ',zqm(igcm_h2o_vap)
	342	print*,'jadrs=',jadrs
	343
	344	call abort
	345	endif
	346	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	347
	348
[135]	349	ENDDO
[253]	350
[3236]	351	pdhadj(:,:)=(zh2(:,:)-zh(:,:))/ptimestep
	352	pduadj(:,:)=(zu2(:,:)-zu(:,:))/ptimestep
	353	pdvadj(:,:)=(zv2(:,:)-zv(:,:))/ptimestep
[135]	354
	355	if(tracer) then
[3236]	356	pdqadj(:,:,:)=(zq2(:,:,:)-zq(:,:,:))/ptimestep
[135]	357	end if
	358
[3322]	359	end subroutine convadj
[253]	360
[3322]	361	end module convadj_mod

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

Download in other formats:

Original Format