Context Navigation

convadj.F @ 2616

Last change on this file since 2616 was 2586, checked in by romain.vande, 3 years ago
LMDZ_MARS Implementation of Open_MP in the physic. Works with radiative transfert
File size: 11.7 KB

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

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

Download in other formats:

Original Format