Context Navigation

moistadj.F90 @ 524

Last change on this file since 524 was 253, checked in by emillour, 14 years ago

Generic GCM

Massive update to version 0.7

EM+RW

File size: 12.3 KB

Rev	Line
[253]	1	subroutine moistadj(t, pq, pplev, pplay, dtmana, dqmana, ptimestep, rneb)
[135]	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
[253]	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	new_qb(k)=0.0
[135]	120	ENDDO
	121
	122	! Calculate v_cptt
	123	DO k = 1, nlayermx
	124	DO i = 1, ngridmx
	125	v_cptt(i,k) = RCPD * local_t(i,k)
	126	ENDDO
	127	ENDDO
	128
	129	DO k = 1, nlayermx
	130	DO i = 1, ngridmx
	131	v_cptt(i,k) = RCPD * local_t(i,k)
	132	v_t = local_t(i,k)
	133	v_p = pplay(i,k)
	134
[253]	135	call watersat(v_t,v_p,v_qs(i,k))
[135]	136	call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
	137	ENDDO
	138	ENDDO
	139
	140	! TEST: RH DIAGNOSTIC
	141	! DO k = 1, nlayermx
	142	! DO i = 1, ngridmx
	143	! v_t = local_t(i,k)
[253]	144	! IF (v_t.LT.To) THEN
[135]	145	! print*,'RHs=',q(i,k) / v_qs(i,k)
	146	! ELSE
	147	! print*,'RHl=',q(i,k) / v_qs(i,k)
	148	! ENDIF
	149	! ENDDO
	150	! ENDDO
	151
	152	! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
	153	DO k = 2, nlayermx
	154	DO i = 1, ngridmx
	155	zdp = pplev(i,k)-pplev(i,k+1)
	156	zdpm = pplev(i,k-1)-pplev(i,k)
	157	! gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
	158	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * &
	159	(v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
	160	+RLVTT /(zdpm+zdp) * &
	161	(v_qs(i,k-1)zdpm + v_qs(i,k)zdp) &
	162	)* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	163	/ (1.0+(v_qsd(i,k-1)*zdpm+ &
	164	v_qsd(i,k)*zdp)/(zdpm+zdp) )
	165	ENDDO
	166	ENDDO
	167
	168	!------------------------------------ modification of unstable profile
	169	DO 9999 i = 1, ngridmx
	170	itest(i) = .FALSE.
	171
	172	! print*,'we in the loop'
	173	! stop
	174
	175	k1 = 0
	176	k2 = 1
	177
	178	810 CONTINUE ! look for k1, the base of the column
	179	k2 = k2 + 1
	180	IF (k2 .GT. nlayermx) GOTO 9999
	181	zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
	182	zsat=(local_q(i,k2-1)-v_qs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) &
	183	+(local_q(i,k2)-v_qs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
	184
	185	IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
	186	k1 = k2 - 1
	187	itest(i) = .TRUE.
	188
	189	820 CONTINUE !! look for k2, the top of the column
	190	IF (k2 .EQ. nlayermx) GOTO 821
	191	k2p = k2 + 1
	192	zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
	193	zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
	194
	195	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821
	196	k2 = k2p
	197	GOTO 820
	198	821 CONTINUE
	199
	200	!------------------------------------------------------ local adjustment
	201	830 CONTINUE ! actual adjustment
	202	v_cptj(k1) = 0.0
	203	zdp = pplev(i,k1)-pplev(i,k1+1)
	204	v_cptjk1 = ( (1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) &
	205	+ RLVTT(local_q(i,k1)-v_qs(i,k1)) ) zdp
	206	v_ssig = zdp * (1.0+v_qsd(i,k1))
	207
	208	k1p = k1 + 1
	209	DO k = k1p, k2
	210	zdp = pplev(i,k)-pplev(i,k+1)
	211	v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
	212	v_cptjk1 = v_cptjk1 + zdp &
	213	* ( (1.0+v_qsd(i, k))*(v_cptt(i,k)+v_cptj(k)) &
	214	+ RLVTT*(local_q(i,k)-v_qs(i,k)) )
	215	v_ssig = v_ssig + zdp *(1.0+v_qsd(i,k))
	216	ENDDO
	217
	218
	219	! this right here is where the adjustment is done???
	220	DO k = k1, k2
	221	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
	222	cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
	223	new_qb(k) = v_qs(i,k) + v_qsd(i,k)*cp_delta_t(k)/RLVTT
	224	local_q(i,k) = new_qb(k)
	225	local_t(i,k) = cp_new_t(k) / RCPD
[253]	226
	227	! print*,'v_qs in loop=',v_qs
	228	! print*,'v_qsd in loop=',v_qsd
	229	! print*,'new_qb in loop=',new_qb
	230	! print*,'cp_delta_t in loop=',cp_delta_t
[135]	231	ENDDO
	232
[253]	233
[135]	234	!--------------------------------------------------- sounding downwards
	235	! -- we refine the prognostic variables in
	236	! -- the layer about to be adjusted
	237
	238	DO k = k1, k2
	239	v_cptt(i,k) = RCPD * local_t(i,k)
	240	v_t = local_t(i,k)
	241	v_p = pplay(i,k)
	242
	243	! IF (v_t.LT.t_coup) THEN
	244	! v_qs(i,k) = qsats(v_t) / v_p
	245	! v_qsd(i,k) = dqsats(v_t,v_qs(i,k))
	246	! ELSE
	247	! v_qs(i,k) = qsatl(v_t) / v_p
	248	! v_qsd(i,k) = dqsatl(v_t,v_qs(i,k))
	249	! ENDIF
	250
[253]	251	call watersat(v_t,v_p,v_qs(i,k))
[135]	252	call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
	253
	254	ENDDO
	255	DO k = 2, nlayermx
	256	zdpm = pplev(i,k-1) - pplev(i,k)
	257	zdp = pplev(i,k) - pplev(i,k+1)
	258	! gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
	259	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * &
	260	(v_cptt(i,k-1)zdpm+v_cptt(i,k)zdp) &
	261	+RLVTT /(zdpm+zdp) * &
	262	(v_qs(i,k-1)zdpm+v_qs(i,k)zdp) &
	263	)* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	264	/ (1.0+(v_qsd(i,k-1)zdpm+v_qsd(i,k)zdp) &
	265	/(zdpm+zdp) )
	266	ENDDO
	267
	268	! Test to see if we've reached the bottom
	269
	270	IF (k1 .EQ. 1) GOTO 841 ! yes we have!
	271	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
	272	zsat=(local_q(i,k1-1)-v_qs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) &
	273	+ (local_q(i,k1)-v_qs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
	274	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
	275
	276	840 CONTINUE
	277	k1 = k1 - 1
	278	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	279	zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1)) &
	280	*(pplev(i,k1-1)-pplev(i,k1))
	281	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
	282	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
	283	GOTO 840
	284	ELSE
	285	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	286	ENDIF
	287	841 CONTINUE
	288
	289	GOTO 810 ! look for other layers higher up
	290
	291	9999 CONTINUE ! loop over all the points
	292
	293	! print*,'k1=',k1
	294	! print*,'k2=',k2
	295
	296	! print*,'local_t=',local_t
	297	! print*,'v_cptt=',v_cptt
	298	! print*,'gamcpdz=',gamcpdz
	299
	300	!-----------------------------------------------------------------------
	301	! Determine the cloud fraction (hypothese: la nebulosite a lieu
	302	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
	303
	304	DO k = 1, nlayermx
	305	DO i = 1, ngridmx
	306	IF (itest(i)) THEN
	307	delta_q(i,k) = local_q(i,k) - q(i,k)
	308	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
	309	ENDIF
	310	ENDDO
	311	ENDDO
	312
	313	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
	314	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
	315	! diminue et d'une maniere proportionnelle a cet diminution):
	316
	317	DO i = 1, ngridmx
	318	IF (itest(i)) THEN
	319	zq1(i) = 0.0
	320	zq2(i) = 0.0
	321	ENDIF
	322	ENDDO
	323	DO k = 1, nlayermx
	324	DO i = 1, ngridmx
	325	IF (itest(i)) THEN
	326	zdp = pplev(i,k)-pplev(i,k+1)
	327	zq1(i) = zq1(i) - delta_q(i,k) * zdp
	328	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
	329	ENDIF
	330	ENDDO
	331	ENDDO
	332	DO k = 1, nlayermx
	333	DO i = 1, ngridmx
	334	IF (itest(i)) THEN
	335	IF (zq2(i).NE.0.0) &
	336	d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
	337	ENDIF
	338	ENDDO
	339	ENDDO
	340
[253]	341	! print*,'local_q BEFORE=',local_q
	342
[135]	343	DO k = 1, nlayermx
	344	DO i = 1, ngridmx
	345	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
	346	ENDDO
	347	ENDDO
	348
	349	DO k = 1, nlayermx
	350	DO i = 1, ngridmx
	351	d_t(i,k) = local_t(i,k) - t(i,k)
	352	d_q(i,k) = local_q(i,k) - q(i,k)
	353	ENDDO
	354	ENDDO
	355
	356	! now subroutine -----> GCM variables
	357	DO k = 1, nlayermx
	358	DO i = 1, ngridmx
	359
	360	dtmana(i,k) = d_t(i,k)/ptimestep
	361	dqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
	362	dqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
	363
	364	ENDDO
	365	ENDDO
	366
	367	! print*,'IN MANABE:'
	368	! print*,'pplev=',pplev
	369	! print*,'t=',t
	370	! print*,'d_t=',d_t
	371	! print*,'d_q=',d_q
[253]	372	! print*,'local_q=',local_q
	373	! print*,'q=',q
	374	! print*,'v_qs(i,k)=',v_qs
	375	! print*,'v_qsd(i,k)=',v_qsd
	376	! print*,'cp_delta_t(k)=',cp_delta_t
	377
[135]	378	! print*,'d_ql=',d_ql
[253]	379	! print*,'delta_q=',delta_q
	380	! print*,'zq1=',zq1
	381	! print*,'zq2=',zq2
	382	!! print*,'i_h2o=',i_h2o
[135]	383	! print*,'i_ice=',i_ice
[253]	384	!
[135]	385	! print*,'IN MANABE:'
	386	! print*,'d_q=',d_q
	387	! print*,'d_ql=',d_ql
	388	! print*,'dtmana=',d_t
	389	! stop
	390	! print*,'gamcpdz at end=',gamcpdz
	391	! stop
	392
	393	! Some conservation diagnostics...
[253]	394	! dEtot=0.0
	395	! dL1tot=0.0
	396	! dL2tot=0.0
	397	! dqtot=0.0
	398	! masse=0.0
	399	! DO k = 1, nlayermx
	400	! DO i = 1, ngridmx
	401	!
	402	! masse = (pplev(i,k) - pplev(i,k+1))/g
	403	!
	404	! dEtot = dEtot + cppd_t(i,k)masse
	405	! dL1tot = dL1tot + RLVTTd_ql(i,k)masse
	406	! dL2tot = dL2tot + RLVTTd_q(i,k)masse ! is this line necessary?
	407	!
	408	! dqtot = dqtot + (d_q(i,k) + d_ql(i,k))*masse
	409	!
	410	! ENDDO
	411	! ENDDO
[135]	412
[253]	413	! print*,'In manabe energy change=',dEtot
	414	! print*,'In manabe condense energy change 1 =',dL1tot
	415	! print*,'In manabe condense energy change 2 =',dL2tot
	416	! print*,'In manabe water change=',dqtot
[135]	417
	418	RETURN
[253]	419	END

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

Download in other formats:

Original Format