Context Navigation

moistadj.F90 @ 1243

Last change on this file since 1243 was 1016, checked in by jleconte, 11 years ago

07/08/2013 == JL

Water cycle in double precision (largescale+moistadj)
Improved wate rayleigh.
First step for rayleigh with variable species. Now, just need to change optcv.
changed some interpolation indices in callcorrk to limit dependency of OLR on the number of layers

File size: 11.2 KB

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

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

Download in other formats:

Original Format