Context Navigation

moistadj.F90 @ 2176

Last change on this file since 2176 was 1384, checked in by emillour, 10 years ago

Generic GCM:

Some code cleanup: turning comcstfi.h into module comcstfi_mod.F90

File size: 11.3 KB

Rev	Line
[1308]	1	subroutine moistadj(ngrid, nlayer, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
[135]	2
[875]	3	use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, RCPV, Psat_water, Lcpdqsat_water
[1308]	4	USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
[1384]	5	use comcstfi_mod, only: r
[135]	6
	7	implicit none
	8
	9
	10	!=====================================================================
	11	!
	12	! Purpose
	13	! -------
	14	! Calculates moist convective adjustment by the method of Manabe.
	15	!
	16	! Authors
	17	! -------
	18	! Adapted from the LMDTERRE code by R. Wordsworth (2010)
	19	! Original author Z. X. Li (1993)
	20	!
	21	!=====================================================================
	22
[1308]	23	INTEGER,INTENT(IN) :: ngrid, nlayer, nq
[135]	24
[1308]	25	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! temperature (K)
	26	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracer (kg/kg)
	27	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
	28	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	29	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
	30	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
	31	REAL,INTENT(OUT) :: pdqmana(ngrid,nlayer,nq) ! tracer tendencies (kg/kg.s-1)
	32	REAL,INTENT(OUT) :: pdtmana(ngrid,nlayer) ! temperature increment(K/s)
	33	REAL,INTENT(OUT) :: rneb(ngrid,nlayer) ! cloud fraction
[135]	34
[728]	35	! local variables
[1308]	36	REAL zt(ngrid,nlayer) ! temperature (K)
	37	REAL zq(ngrid,nlayer) ! humidite specifique (kg/kg)
[135]	38
[1308]	39	REAL d_t(ngrid,nlayer) ! temperature increment
	40	REAL d_q(ngrid,nlayer) ! incrementation pour vapeur d'eau
	41	REAL d_ql(ngrid,nlayer) ! incrementation pour l'eau liquide
[135]	42
	43	! REAL t_coup
	44	! PARAMETER (t_coup=234.0)
	45	REAL seuil_vap
	46	PARAMETER (seuil_vap=1.0E-10)
	47
	48	! Local variables
[728]	49	INTEGER i, k, iq, nn
[773]	50	INTEGER, PARAMETER :: niter = 1
[135]	51	INTEGER k1, k1p, k2, k2p
[787]	52	LOGICAL itest(ngrid)
[1308]	53	REAL delta_q(ngrid, nlayer)
	54	DOUBLE PRECISION :: cp_new_t(nlayer), v_cptt(ngrid,nlayer)
	55	REAL cp_delta_t(nlayer)
	56	DOUBLE PRECISION :: v_cptj(nlayer), v_cptjk1, v_ssig
[1016]	57	REAL v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat
[1308]	58	REAL zqs(ngrid,nlayer), zdqs(ngrid,nlayer),zpsat(ngrid,nlayer),zdlnpsat(ngrid,nlayer)
[787]	59	REAL zq1(ngrid), zq2(ngrid)
[1308]	60	DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayer)
[1016]	61	DOUBLE PRECISION :: zdp, zdpm
[135]	62
	63	REAL zsat ! super-saturation
	64	REAL zflo ! flotabilite
	65
[1308]	66	DOUBLE PRECISION :: local_q(ngrid,nlayer),local_t(ngrid,nlayer)
[135]	67
	68	REAL zdelta, zcor, zcvm5
	69
	70	REAL dEtot, dqtot, masse ! conservation diagnostics
	71	real dL1tot, dL2tot
	72
	73	! Indices of water vapour and water ice tracers
	74	INTEGER,SAVE :: i_h2o=0 ! water vapour
	75	INTEGER,SAVE :: i_ice=0 ! water ice
[1315]	76	!$OMP THREADPRIVATE(i_h2o,i_ice)
[135]	77
[1308]	78	LOGICAL,SAVE :: firstcall=.TRUE.
[1315]	79	!$OMP THREADPRIVATE(firstcall)
[135]	80
	81	IF (firstcall) THEN
	82
	83	i_h2o=igcm_h2o_vap
	84	i_ice=igcm_h2o_ice
	85
	86	write(,) "rain: i_ice=",i_ice
	87	write(,) " i_h2o=",i_h2o
	88
	89	firstcall = .FALSE.
	90	ENDIF
	91
	92	! GCM -----> subroutine variables
[1308]	93	zq(1:ngrid,1:nlayer) = pq(1:ngrid,1:nlayer,i_h2o)+ pdq(1:ngrid,1:nlayer,i_h2o)*ptimestep
	94	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer)
	95	pdqmana(1:ngrid,1:nlayer,1:nq)=0.0
[728]	96
[1308]	97	DO k = 1, nlayer
[787]	98	DO i = 1, ngrid
[728]	99	if(zq(i,k).lt.0.)then
	100	zq(i,k)=0.0
[135]	101	endif
[728]	102	ENDDO
[135]	103	ENDDO
[728]	104
[1308]	105	local_q(1:ngrid,1:nlayer) = zq(1:ngrid,1:nlayer)
	106	local_t(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
	107	rneb(1:ngrid,1:nlayer) = 0.0
	108	d_ql(1:ngrid,1:nlayer) = 0.0
	109	d_t(1:ngrid,1:nlayer) = 0.0
	110	d_q(1:ngrid,1:nlayer) = 0.0
[135]	111
	112	! Calculate v_cptt
[1308]	113	DO k = 1, nlayer
[787]	114	DO i = 1, ngrid
[135]	115	v_cptt(i,k) = RCPD * local_t(i,k)
[728]	116	v_t = MAX(local_t(i,k),15.)
[135]	117	v_p = pplay(i,k)
	118
[875]	119	call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
	120	call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
[135]	121	ENDDO
	122	ENDDO
	123
	124	! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
[1308]	125	DO k = 2, nlayer
[787]	126	DO i = 1, ngrid
[135]	127	zdp = pplev(i,k)-pplev(i,k+1)
	128	zdpm = pplev(i,k-1)-pplev(i,k)
[875]	129	! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
	130	! + RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
	131	! * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	132	! / (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
	133	! general case where water is not a trace gas (JL13)
	134	v_zqs = (zqs(i,k-1)zdpm + zqs(i,k)zdp)/(zdpm+zdp)
	135	v_cptt2 = (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp)/(zdpm+zdp)
[1016]	136	v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))zdpm + (1.-zpsat(i,k)/pplay(i,k))zdp)/(zdpm+zdp)
[875]	137	v_dlnpsat = (zdlnpsat(i,k-1)zdpm + zdlnpsat(i,k)zdp)/(zdpm+zdp)
[1016]	138	gamcpdz(i,k) = ( (R/RCPDv_cptt2(1.- v_zqs) + RLVTTv_zqs) (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) &
	139	/ (((1.- v_zqs) + v_zqs * RCPV/RCPD)v_pratio + v_zqs v_dlnpsat)
[135]	140	ENDDO
	141	ENDDO
	142
	143	!------------------------------------ modification of unstable profile
[787]	144	DO 9999 i = 1, ngrid
	145
[135]	146	itest(i) = .FALSE.
	147
	148	! print*,'we in the loop'
	149	! stop
	150
	151	k1 = 0
	152	k2 = 1
	153
	154	810 CONTINUE ! look for k1, the base of the column
	155	k2 = k2 + 1
[1308]	156	IF (k2 .GT. nlayer) GOTO 9999
[135]	157	zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
[728]	158	zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) &
	159	+(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
[135]	160
	161	IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
	162	k1 = k2 - 1
	163	itest(i) = .TRUE.
	164
	165	820 CONTINUE !! look for k2, the top of the column
[1308]	166	IF (k2 .EQ. nlayer) GOTO 821
[135]	167	k2p = k2 + 1
[728]	168	zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
[135]	169	zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
	170
	171	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821
	172	k2 = k2p
	173	GOTO 820
	174	821 CONTINUE
	175
	176	!------------------------------------------------------ local adjustment
	177	830 CONTINUE ! actual adjustment
[728]	178	Do nn=1,niter
[135]	179	v_cptj(k1) = 0.0
	180	zdp = pplev(i,k1)-pplev(i,k1+1)
[728]	181	v_cptjk1 = ( (1.0+zdqs(i,k1))(v_cptt(i,k1)+v_cptj(k1)) + RLVTT(local_q(i,k1)-zqs(i,k1)) ) * zdp
	182	v_ssig = zdp * (1.0+zdqs(i,k1))
[135]	183
	184	k1p = k1 + 1
	185	DO k = k1p, k2
	186	zdp = pplev(i,k)-pplev(i,k+1)
	187	v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
[728]	188	v_cptjk1 = v_cptjk1 + zdp * ( (1.0+zdqs(i, k))(v_cptt(i,k)+v_cptj(k)) + RLVTT(local_q(i,k)-zqs(i,k)) )
	189	v_ssig = v_ssig + zdp *(1.0+zdqs(i,k))
[135]	190	ENDDO
	191
	192
	193	! this right here is where the adjustment is done???
	194	DO k = k1, k2
	195	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
	196	cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
[728]	197	v_cptt(i,k)=cp_new_t(k)
	198	local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT
[135]	199	local_t(i,k) = cp_new_t(k) / RCPD
[253]	200
[728]	201	v_t = MAX(local_t(i,k),15.)
	202	v_p = pplay(i,k)
	203
[875]	204	call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
	205	call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
[728]	206
[135]	207	ENDDO
[728]	208	Enddo ! nn=1,niter
[135]	209
[253]	210
[135]	211	!--------------------------------------------------- sounding downwards
	212	! -- we refine the prognostic variables in
	213	! -- the layer about to be adjusted
	214
[728]	215	! DO k = k1, k2
	216	! v_cptt(i,k) = RCPD * local_t(i,k)
	217	! v_t = local_t(i,k)
	218	! v_p = pplay(i,k)
	219	!
	220	! call Psat_water(v_t,v_p,zpsat,zqs(i,k))
	221	! call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
	222	! ENDDO
[135]	223
[1308]	224	DO k = 2, nlayer
[135]	225	zdpm = pplev(i,k-1) - pplev(i,k)
	226	zdp = pplev(i,k) - pplev(i,k+1)
[875]	227	! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
	228	! + RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
	229	! * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	230	! / (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
	231	! general case where water is not a trace gas (JL13)
	232	v_zqs = (zqs(i,k-1)zdpm + zqs(i,k)zdp)/(zdpm+zdp)
	233	v_cptt2 = (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp)/(zdpm+zdp)
[1016]	234	v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))zdpm + (1.-zpsat(i,k)/pplay(i,k))zdp)/(zdpm+zdp)
[875]	235	v_dlnpsat = (zdlnpsat(i,k-1)zdpm + zdlnpsat(i,k)zdp)/(zdpm+zdp)
[1016]	236	gamcpdz(i,k) = ( (R/RCPDv_cptt2(1.- v_zqs) + RLVTTv_zqs) (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) &
	237	/ (((1.- v_zqs) + v_zqs * RCPV/RCPD)v_pratio + v_zqs v_dlnpsat)
[135]	238	ENDDO
	239
	240	! Test to see if we've reached the bottom
	241
	242	IF (k1 .EQ. 1) GOTO 841 ! yes we have!
	243	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
[728]	244	zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) &
	245	+ (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
[135]	246	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
	247
	248	840 CONTINUE
	249	k1 = k1 - 1
	250	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
[728]	251	zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1)) &
[135]	252	*(pplev(i,k1-1)-pplev(i,k1))
	253	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
	254	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
	255	GOTO 840
	256	ELSE
	257	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	258	ENDIF
	259	841 CONTINUE
	260
	261	GOTO 810 ! look for other layers higher up
	262
	263	9999 CONTINUE ! loop over all the points
	264
	265	!-----------------------------------------------------------------------
	266	! Determine the cloud fraction (hypothese: la nebulosite a lieu
	267	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
	268
[1308]	269	DO k = 1, nlayer
[787]	270	DO i = 1, ngrid
[135]	271	IF (itest(i)) THEN
[728]	272	delta_q(i,k) = local_q(i,k) - zq(i,k)
[135]	273	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
	274	ENDIF
	275	ENDDO
	276	ENDDO
	277
	278	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
	279	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
	280	! diminue et d'une maniere proportionnelle a cet diminution):
	281
[787]	282	DO i = 1, ngrid
[135]	283	IF (itest(i)) THEN
	284	zq1(i) = 0.0
	285	zq2(i) = 0.0
	286	ENDIF
	287	ENDDO
[1308]	288	DO k = 1, nlayer
[787]	289	DO i = 1, ngrid
[135]	290	IF (itest(i)) THEN
	291	zdp = pplev(i,k)-pplev(i,k+1)
	292	zq1(i) = zq1(i) - delta_q(i,k) * zdp
	293	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
	294	ENDIF
	295	ENDDO
	296	ENDDO
[1308]	297	DO k = 1, nlayer
[787]	298	DO i = 1, ngrid
[135]	299	IF (itest(i)) THEN
[728]	300	IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
[135]	301	ENDIF
	302	ENDDO
	303	ENDDO
	304
[1308]	305	DO k = 1, nlayer
[787]	306	DO i = 1, ngrid
[135]	307	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
	308	ENDDO
	309	ENDDO
	310
[1308]	311	DO k = 1, nlayer
[787]	312	DO i = 1, ngrid
[728]	313	d_t(i,k) = local_t(i,k) - zt(i,k)
	314	d_q(i,k) = local_q(i,k) - zq(i,k)
[135]	315	ENDDO
	316	ENDDO
	317
	318	! now subroutine -----> GCM variables
[1308]	319	DO k = 1, nlayer
[787]	320	DO i = 1, ngrid
[135]	321
[728]	322	pdtmana(i,k) = d_t(i,k)/ptimestep
	323	pdqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
	324	pdqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
[135]	325
	326	ENDDO
	327	ENDDO
	328
[253]	329
	330	END

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

Download in other formats:

Original Format