Context Navigation

moistadj.F90 @ 164

Last change on this file since 164 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 11.8 KB

Rev	Line
[135]	1	subroutine moistadj(t, pq, pplev, pplay, dtmana, dqmana, ptimestep, rneb)
	2
	3	use watercommon_h, only: To, RLVTT, RCPD
	4
	5	implicit none
	6
	7
	8	!=====================================================================
	9	!
	10	! Purpose
	11	! -------
	12	! Calculates moist convective adjustment by the method of Manabe.
	13	!
	14	! Authors
	15	! -------
	16	! Adapted from the LMDTERRE code by R. Wordsworth (2010)
	17	! Original author Z. X. Li (1993)
	18	!
	19	!=====================================================================
	20
	21	#include "dimensions.h"
	22	#include "dimphys.h"
	23	#include "tracer.h"
	24	#include "comcstfi.h"
	25
	26	! Pre-arguments (for universal model)
	27	real pq(ngridmx,nlayermx,nqmx) ! tracer (kg/kg)
	28	REAL pdq(ngridmx,nlayermx,nqmx)
	29
	30	real dqmana(ngridmx,nlayermx,nqmx) ! tendency of tracers (kg/kg.s-1)
	31	REAL dtmana(ngridmx,nlayermx) ! temperature increment
	32
	33	! Arguments
	34	REAL t(ngridmx,nlayermx) ! temperature (K)
	35	REAL q(ngridmx,nlayermx) ! humidite specifique (kg/kg)
	36	REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche (Pa)
	37	REAL pplay(ngridmx,nlayermx) ! pression au milieu de couche (Pa)
	38
	39	REAL d_t(ngridmx,nlayermx) ! temperature increment
	40	REAL d_q(ngridmx,nlayermx) ! incrementation pour vapeur d'eau
	41	REAL d_ql(ngridmx,nlayermx) ! incrementation pour l'eau liquide
	42	REAL rneb(ngridmx,nlayermx) ! cloud fraction
	43	REAL ptimestep
	44
	45	! REAL t_coup
	46	! PARAMETER (t_coup=234.0)
	47	REAL seuil_vap
	48	PARAMETER (seuil_vap=1.0E-10)
	49
	50	! Local variables
	51	INTEGER i, k, iq
	52	INTEGER k1, k1p, k2, k2p
	53	LOGICAL itest(ngridmx)
	54	REAL delta_q(ngridmx, nlayermx)
	55	REAL cp_new_t(nlayermx)
	56	REAL cp_delta_t(nlayermx)
	57	REAL new_qb(nlayermx)
	58	REAL v_cptj(nlayermx), v_cptjk1, v_ssig
	59	REAL v_cptt(ngridmx,nlayermx), v_p, v_t
	60	REAL v_qs(ngridmx,nlayermx), v_qsd(ngridmx,nlayermx)
	61	REAL zq1(ngridmx), zq2(ngridmx)
	62	REAL gamcpdz(ngridmx,2:nlayermx)
	63	REAL zdp, zdpm
	64
	65	REAL zsat ! super-saturation
	66	REAL zflo ! flotabilite
	67
	68	REAL local_q(ngridmx,nlayermx),local_t(ngridmx,nlayermx)
	69
	70	REAL zdelta, zcor, zcvm5
	71
	72	REAL dEtot, dqtot, masse ! conservation diagnostics
	73	real dL1tot, dL2tot
	74
	75	! Indices of water vapour and water ice tracers
	76	INTEGER,SAVE :: i_h2o=0 ! water vapour
	77	INTEGER,SAVE :: i_ice=0 ! water ice
	78
	79	LOGICAL firstcall
	80	SAVE firstcall
	81
	82	DATA firstcall /.TRUE./
	83
	84	IF (firstcall) THEN
	85
	86	i_h2o=igcm_h2o_vap
	87	i_ice=igcm_h2o_ice
	88
	89	write(,) "rain: i_ice=",i_ice
	90	write(,) " i_h2o=",i_h2o
	91
	92	firstcall = .FALSE.
	93	ENDIF
	94
	95	! GCM -----> subroutine variables
	96	DO k = 1, nlayermx
	97	DO i = 1, ngridmx
	98
	99	q(i,k) = pq(i,k,i_h2o)
	100
	101	if(q(i,k).lt.0.)then
	102	q(i,k)=0.0
	103	endif
	104	DO iq = 1, nqmx
	105	dqmana(i,k,iq)=0.0
	106	ENDDO
	107	ENDDO
	108	ENDDO
	109
	110	DO k = 1, nlayermx
	111	DO i = 1, ngridmx
	112	local_q(i,k) = q(i,k)
	113	local_t(i,k) = t(i,k)
	114	rneb(i,k) = 0.0
	115	d_ql(i,k) = 0.0
	116	d_t(i,k) = 0.0
	117	d_q(i,k) = 0.0
	118	ENDDO
	119	ENDDO
	120
	121	! Calculate v_cptt
	122	DO k = 1, nlayermx
	123	DO i = 1, ngridmx
	124	v_cptt(i,k) = RCPD * local_t(i,k)
	125	ENDDO
	126	ENDDO
	127
	128	DO k = 1, nlayermx
	129	DO i = 1, ngridmx
	130	v_cptt(i,k) = RCPD * local_t(i,k)
	131	v_t = local_t(i,k)
	132	v_p = pplay(i,k)
	133
	134	call watersat_2(v_t,v_p,v_qs(i,k))
	135	call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
	136	ENDDO
	137	ENDDO
	138
	139	! TEST: RH DIAGNOSTIC
	140	! DO k = 1, nlayermx
	141	! DO i = 1, ngridmx
	142	! v_t = local_t(i,k)
	143	! IF (v_t.LT.t_coup) THEN
	144	! print*,'RHs=',q(i,k) / v_qs(i,k)
	145	! ELSE
	146	! print*,'RHl=',q(i,k) / v_qs(i,k)
	147	! ENDIF
	148	! ENDDO
	149	! ENDDO
	150
	151	! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
	152	DO k = 2, nlayermx
	153	DO i = 1, ngridmx
	154	zdp = pplev(i,k)-pplev(i,k+1)
	155	zdpm = pplev(i,k-1)-pplev(i,k)
	156	! gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
	157	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * &
	158	(v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
	159	+RLVTT /(zdpm+zdp) * &
	160	(v_qs(i,k-1)zdpm + v_qs(i,k)zdp) &
	161	)* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	162	/ (1.0+(v_qsd(i,k-1)*zdpm+ &
	163	v_qsd(i,k)*zdp)/(zdpm+zdp) )
	164	ENDDO
	165	ENDDO
	166
	167	!------------------------------------ modification of unstable profile
	168	DO 9999 i = 1, ngridmx
	169	itest(i) = .FALSE.
	170
	171	! print*,'we in the loop'
	172	! stop
	173
	174	k1 = 0
	175	k2 = 1
	176
	177	810 CONTINUE ! look for k1, the base of the column
	178	k2 = k2 + 1
	179	IF (k2 .GT. nlayermx) GOTO 9999
	180	zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
	181	zsat=(local_q(i,k2-1)-v_qs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) &
	182	+(local_q(i,k2)-v_qs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
	183
	184	IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
	185	k1 = k2 - 1
	186	itest(i) = .TRUE.
	187
	188	820 CONTINUE !! look for k2, the top of the column
	189	IF (k2 .EQ. nlayermx) GOTO 821
	190	k2p = k2 + 1
	191	zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
	192	zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
	193
	194	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821
	195	k2 = k2p
	196	GOTO 820
	197	821 CONTINUE
	198
	199	!------------------------------------------------------ local adjustment
	200	830 CONTINUE ! actual adjustment
	201	v_cptj(k1) = 0.0
	202	zdp = pplev(i,k1)-pplev(i,k1+1)
	203	v_cptjk1 = ( (1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) &
	204	+ RLVTT(local_q(i,k1)-v_qs(i,k1)) ) zdp
	205	v_ssig = zdp * (1.0+v_qsd(i,k1))
	206
	207	k1p = k1 + 1
	208	DO k = k1p, k2
	209	zdp = pplev(i,k)-pplev(i,k+1)
	210	v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
	211	v_cptjk1 = v_cptjk1 + zdp &
	212	* ( (1.0+v_qsd(i, k))*(v_cptt(i,k)+v_cptj(k)) &
	213	+ RLVTT*(local_q(i,k)-v_qs(i,k)) )
	214	v_ssig = v_ssig + zdp *(1.0+v_qsd(i,k))
	215	ENDDO
	216
	217
	218	! this right here is where the adjustment is done???
	219	DO k = k1, k2
	220	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
	221	cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
	222	new_qb(k) = v_qs(i,k) + v_qsd(i,k)*cp_delta_t(k)/RLVTT
	223	local_q(i,k) = new_qb(k)
	224	local_t(i,k) = cp_new_t(k) / RCPD
	225	ENDDO
	226
	227	!--------------------------------------------------- sounding downwards
	228	! -- we refine the prognostic variables in
	229	! -- the layer about to be adjusted
	230
	231	DO k = k1, k2
	232	v_cptt(i,k) = RCPD * local_t(i,k)
	233	v_t = local_t(i,k)
	234	v_p = pplay(i,k)
	235
	236	! IF (v_t.LT.t_coup) THEN
	237	! v_qs(i,k) = qsats(v_t) / v_p
	238	! v_qsd(i,k) = dqsats(v_t,v_qs(i,k))
	239	! ELSE
	240	! v_qs(i,k) = qsatl(v_t) / v_p
	241	! v_qsd(i,k) = dqsatl(v_t,v_qs(i,k))
	242	! ENDIF
	243
	244	call watersat_2(v_t,v_p,v_qs(i,k))
	245	call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
	246
	247	ENDDO
	248	DO k = 2, nlayermx
	249	zdpm = pplev(i,k-1) - pplev(i,k)
	250	zdp = pplev(i,k) - pplev(i,k+1)
	251	! gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
	252	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * &
	253	(v_cptt(i,k-1)zdpm+v_cptt(i,k)zdp) &
	254	+RLVTT /(zdpm+zdp) * &
	255	(v_qs(i,k-1)zdpm+v_qs(i,k)zdp) &
	256	)* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	257	/ (1.0+(v_qsd(i,k-1)zdpm+v_qsd(i,k)zdp) &
	258	/(zdpm+zdp) )
	259	ENDDO
	260
	261	! Test to see if we've reached the bottom
	262
	263	IF (k1 .EQ. 1) GOTO 841 ! yes we have!
	264	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
	265	zsat=(local_q(i,k1-1)-v_qs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) &
	266	+ (local_q(i,k1)-v_qs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
	267	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
	268
	269	840 CONTINUE
	270	k1 = k1 - 1
	271	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	272	zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1)) &
	273	*(pplev(i,k1-1)-pplev(i,k1))
	274	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
	275	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
	276	GOTO 840
	277	ELSE
	278	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	279	ENDIF
	280	841 CONTINUE
	281
	282	GOTO 810 ! look for other layers higher up
	283
	284	9999 CONTINUE ! loop over all the points
	285
	286	! print*,'k1=',k1
	287	! print*,'k2=',k2
	288
	289	! print*,'local_t=',local_t
	290	! print*,'v_cptt=',v_cptt
	291	! print*,'gamcpdz=',gamcpdz
	292
	293	!-----------------------------------------------------------------------
	294	! Determine the cloud fraction (hypothese: la nebulosite a lieu
	295	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
	296
	297	DO k = 1, nlayermx
	298	DO i = 1, ngridmx
	299	IF (itest(i)) THEN
	300	delta_q(i,k) = local_q(i,k) - q(i,k)
	301	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
	302	ENDIF
	303	ENDDO
	304	ENDDO
	305
	306	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
	307	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
	308	! diminue et d'une maniere proportionnelle a cet diminution):
	309
	310	DO i = 1, ngridmx
	311	IF (itest(i)) THEN
	312	zq1(i) = 0.0
	313	zq2(i) = 0.0
	314	ENDIF
	315	ENDDO
	316	DO k = 1, nlayermx
	317	DO i = 1, ngridmx
	318	IF (itest(i)) THEN
	319	zdp = pplev(i,k)-pplev(i,k+1)
	320	zq1(i) = zq1(i) - delta_q(i,k) * zdp
	321	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
	322	ENDIF
	323	ENDDO
	324	ENDDO
	325	DO k = 1, nlayermx
	326	DO i = 1, ngridmx
	327	IF (itest(i)) THEN
	328	IF (zq2(i).NE.0.0) &
	329	d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
	330	ENDIF
	331	ENDDO
	332	ENDDO
	333
	334	DO k = 1, nlayermx
	335	DO i = 1, ngridmx
	336	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
	337	ENDDO
	338	ENDDO
	339
	340	DO k = 1, nlayermx
	341	DO i = 1, ngridmx
	342	d_t(i,k) = local_t(i,k) - t(i,k)
	343	d_q(i,k) = local_q(i,k) - q(i,k)
	344	ENDDO
	345	ENDDO
	346
	347	! now subroutine -----> GCM variables
	348	DO k = 1, nlayermx
	349	DO i = 1, ngridmx
	350
	351	dtmana(i,k) = d_t(i,k)/ptimestep
	352	dqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
	353	dqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
	354
	355	ENDDO
	356	ENDDO
	357
	358	! print*,'IN MANABE:'
	359	! print*,'pplev=',pplev
	360	! print*,'t=',t
	361	! print*,'d_t=',d_t
	362	! print*,'d_q=',d_q
	363	! print*,'d_ql=',d_ql
	364	! print*,'i_h2o=',i_h2o
	365	! print*,'i_ice=',i_ice
	366
	367	! print*,'IN MANABE:'
	368	! print*,'d_q=',d_q
	369	! print*,'d_ql=',d_ql
	370	! print*,'dtmana=',d_t
	371	! stop
	372	! print*,'gamcpdz at end=',gamcpdz
	373	! stop
	374
	375	! Some conservation diagnostics...
	376	dEtot=0.0
	377	dL1tot=0.0
	378	dL2tot=0.0
	379	dqtot=0.0
	380	masse=0.0
	381	DO k = 1, nlayermx
	382	DO i = 1, ngridmx
	383
	384	masse = (pplev(i,k) - pplev(i,k+1))/g
	385
	386	dEtot = dEtot + cppd_t(i,k)masse
	387	dL1tot = dL1tot + RLVTTd_ql(i,k)masse
	388	dL2tot = dL2tot + RLVTTd_q(i,k)masse ! is this line necessary?
	389
	390	dqtot = dqtot + (d_q(i,k) + d_ql(i,k))*masse
	391
	392	ENDDO
	393	ENDDO
	394
	395	print*,'In manabe energy change=',dEtot
	396	print*,'In manabe condense energy change 1 =',dL1tot
	397	print*,'In manabe condense energy change 2 =',dL2tot
	398	print*,'In manabe water change=',dqtot
	399
	400	RETURN
	401	END

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

Download in other formats:

Original Format