Context Navigation

moistadj_generic.F90 @ 3105

Last change on this file since 3105 was 3104, checked in by emillour, 20 months ago
new routine moistadj_generic
File size: 14.0 KB

Rev	Line
[3104]	1	subroutine moistadj_generic(ngrid, nlayer, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
	2
	3	!use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, RCPV, Psat_water, Lcpdqsat_water
	4	!USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
	5	use generic_cloud_common_h
	6	use generic_tracer_index_mod, only: generic_tracer_index
	7	use tracer_h
	8	use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity
	9	use comcstfi_mod, only: r, cpp, mugaz
	10
	11	implicit none
	12
	13
	14	!=====================================================================
	15	!
	16	! Purpose
	17	! -------
	18	! Calculates moist convective adjustment by the method of Manabe.
	19	!
	20	! Authors
	21	! -------
	22	! Adapted from the moistadj.F90 routine
	23	! for generic tracers (condensable species)
	24	! by Noe CLEMENT (2023)
	25	!
	26	!=====================================================================
	27
	28	INTEGER,INTENT(IN) :: ngrid, nlayer, nq
	29
	30	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! temperature (K)
	31	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracer (kg/kg)
	32	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
	33	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	34	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
	35	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
	36	REAL,INTENT(OUT) :: pdqmana(ngrid,nlayer,nq) ! tracer tendencies (kg/kg.s-1)
	37	REAL,INTENT(OUT) :: pdtmana(ngrid,nlayer) ! temperature increment(K/s)
	38	REAL,INTENT(OUT) :: rneb(ngrid,nlayer) ! cloud fraction
	39
	40	! Options :
	41	!real, save :: metallicity ! metallicity of planet
	42	REAL,SAVE :: metallicity = 0.0
	43	!$OMP THREADPRIVATE(metallicity)
	44
	45	! local variables
	46	REAL zt(ngrid,nlayer) ! temperature (K)
	47	REAL zq(ngrid,nlayer) ! humidite specifique (kg/kg)
	48
	49	REAL d_t(ngrid,nlayer) ! temperature increment
	50	REAL d_q(ngrid,nlayer) ! incrementation pour vapeur d'eau
	51	REAL d_ql(ngrid,nlayer) ! incrementation pour l'eau liquide
	52
	53	! REAL t_coup
	54	! PARAMETER (t_coup=234.0)
	55	REAL seuil_vap
	56	PARAMETER (seuil_vap=1.0E-10)
	57
	58	! Local variables
	59
	60	logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
	61	integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice of generic_tracer
	62
	63	INTEGER i, k, iq, nn
	64	INTEGER, PARAMETER :: niter = 1
	65	INTEGER k1, k1p, k2, k2p
	66	LOGICAL itest(ngrid)
	67	REAL delta_q(ngrid, nlayer)
	68	DOUBLE PRECISION :: cp_new_t(nlayer), v_cptt(ngrid,nlayer)
	69	REAL cp_delta_t(nlayer)
	70	DOUBLE PRECISION :: v_cptj(nlayer), v_cptjk1, v_ssig
	71	REAL v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat
	72	REAL zqs(ngrid,nlayer), zdqs(ngrid,nlayer),zpsat(ngrid,nlayer),zdlnpsat(ngrid,nlayer)
	73	REAL zq1(ngrid), zq2(ngrid)
	74	DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayer)
	75	DOUBLE PRECISION :: zdp, zdpm
	76
	77	real q_cri(ngrid,nlayer)
	78
	79	REAL zsat ! super-saturation
	80	REAL zflo ! flotabilite
	81
	82	DOUBLE PRECISION :: local_q(ngrid,nlayer),local_t(ngrid,nlayer)
	83
	84	REAL zdelta, zcor, zcvm5
	85
	86	REAL dEtot, dqtot, masse ! conservation diagnostics
	87	real dL1tot, dL2tot
	88
	89	! Indices of generic vapour and ice tracers
	90	real,save :: RCPD=0.0
	91	INTEGER,SAVE :: i_vap_generic=0 ! Generic Condensable Species vapour
	92	INTEGER,SAVE :: i_ice_generic=0 ! Generic Condensable Species ice
	93	!$OMP THREADPRIVATE(i_vap_generic,i_ice_generic,RCPD)
	94
	95	LOGICAL,SAVE :: firstcall=.TRUE.
	96	!$OMP THREADPRIVATE(firstcall)
	97
	98	IF (firstcall) THEN
	99
	100	RCPD = cpp
	101	! RCPV = 1.88e3 ! specific heat capacity of water vapor at 350K ! MUST BE CHANGED !!!
	102	! 1.709e3 for methane - should be added in table_tracers_condensable (datagcm)
	103
	104	write(,) "value for metallicity? "
	105	metallicity=0.0 ! default value
	106	call getin_p("metallicity",metallicity)
	107	write(,) " metallicity = ",metallicity
	108
	109	do iq=1, nq
	110	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
	111	if(call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
	112	write(,) "moistadjustment : taking back the values you have set in 'table_tracers_condensable' for thermodynamics. If you have forgotten one, it will crash."
	113	m = constants_mass(iq)
	114	delta_vapH = constants_delta_vapH(iq)
	115	epsi_generic = constants_epsi_generic(iq)
	116	RLVTT_generic = constants_RLVTT_generic(iq)
	117	RCPV_generic = constants_RCPV_generic(iq)
	118	endif
	119	enddo
	120
	121	i_vap_generic=igcm_generic_vap
	122	i_ice_generic=igcm_generic_ice
	123
	124	write(,) "rain: i_ice_generic=",i_ice_generic
	125	write(,) " i_vap_generic=",i_vap_generic
	126
	127	firstcall = .FALSE.
	128	ENDIF
	129
	130	! GCM -----> subroutine variables
	131	zq(1:ngrid,1:nlayer) = pq(1:ngrid,1:nlayer,i_vap_generic)+ pdq(1:ngrid,1:nlayer,i_vap_generic)*ptimestep
	132	zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer)
	133	pdqmana(1:ngrid,1:nlayer,1:nq)=0.0
	134
	135	DO k = 1, nlayer
	136	DO i = 1, ngrid
	137	if(zq(i,k).lt.0.)then
	138	zq(i,k)=0.0
	139	endif
	140	ENDDO
	141	ENDDO
	142
	143	local_q(1:ngrid,1:nlayer) = zq(1:ngrid,1:nlayer)
	144	local_t(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
	145	rneb(1:ngrid,1:nlayer) = 0.0
	146	d_ql(1:ngrid,1:nlayer) = 0.0
	147	d_t(1:ngrid,1:nlayer) = 0.0
	148	d_q(1:ngrid,1:nlayer) = 0.0
	149
	150	! Calculate v_cptt
	151	DO k = 1, nlayer
	152	DO i = 1, ngrid
	153	v_cptt(i,k) = RCPD * local_t(i,k)
	154	v_t = MAX(local_t(i,k),15.)
	155	v_p = pplay(i,k)
	156
	157	call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k))
	158	call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
	159
	160	!call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
	161	!call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
	162	ENDDO
	163	ENDDO
	164
	165	! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
	166	DO k = 2, nlayer
	167	DO i = 1, ngrid
	168	zdp = pplev(i,k)-pplev(i,k+1)
	169	zdpm = pplev(i,k-1)-pplev(i,k)
	170	! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
	171	! + RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
	172	!* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	173	! / (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
	174	! general case where water is not a trace gas (JL13)
	175	v_zqs = (zqs(i,k-1)zdpm + zqs(i,k)zdp)/(zdpm+zdp)
	176	v_cptt2 = (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp)/(zdpm+zdp)
	177	v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))zdpm + (1.-zpsat(i,k)/pplay(i,k))zdp)/(zdpm+zdp)
	178	v_dlnpsat = (zdlnpsat(i,k-1)zdpm + zdlnpsat(i,k)zdp)/(zdpm+zdp)
	179	gamcpdz(i,k) = ( (R/RCPDv_cptt2(1.- v_zqs) + RLVTT_genericv_zqs) (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) &
	180	/ (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)v_pratio + v_zqs v_dlnpsat)
	181	ENDDO
	182	ENDDO
	183
	184	! calculates moist convection inhibition criterion
	185	DO k = 1, nlayer
	186	DO i = 1, ngrid
	187	q_cri(i,k) = ( 1 / (1 - 1/epsi_generic)) * r * mugaz/1000. / delta_vapH * zt(i,k)
	188	ENDDO
	189	ENDDO
	190
	191	write(,) "q_cri at 80K : ", ( 1 / (1 - 1/epsi_generic)) * (r * mugaz/1000.) / delta_vapH * 80
	192
	193	!------------------------------------ modification of unstable profile
	194	DO 9999 i = 1, ngrid
	195
	196	itest(i) = .FALSE.
	197
	198	! print*,'we in the loop'
	199	! stop
	200
	201	k1 = 0
	202	k2 = 1
	203
	204	810 CONTINUE ! look for k1, the base of the column
	205	k2 = k2 + 1
	206	IF (k2 .GT. nlayer) GOTO 9999
	207	zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
	208	zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) &
	209	+(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
	210
	211	IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 .OR. local_q(i,k2-1).GT.q_cri(i,k2-1)) GOTO 810
	212	k1 = k2 - 1
	213	itest(i) = .TRUE.
	214
	215	820 CONTINUE !! look for k2, the top of the column
	216	IF (k2 .EQ. nlayer) GOTO 821
	217	k2p = k2 + 1
	218	zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
	219	zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
	220
	221	IF (zflo.LE.0.0 .OR. zsat.LE.0.0 .OR. local_q(i,k2p).GT.q_cri(i,k2p)) GOTO 821
	222	k2 = k2p
	223	GOTO 820
	224	821 CONTINUE
	225
	226	!------------------------------------------------------ local adjustment
	227	830 CONTINUE ! actual adjustment
	228	Do nn=1,niter
	229	v_cptj(k1) = 0.0
	230	zdp = pplev(i,k1)-pplev(i,k1+1)
	231	v_cptjk1 = ( (1.0+zdqs(i,k1))(v_cptt(i,k1)+v_cptj(k1)) + RLVTT_generic(local_q(i,k1)-zqs(i,k1)) ) * zdp
	232	v_ssig = zdp * (1.0+zdqs(i,k1))
	233
	234	k1p = k1 + 1
	235	DO k = k1p, k2
	236	zdp = pplev(i,k)-pplev(i,k+1)
	237	v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
	238	v_cptjk1 = v_cptjk1 + zdp * ( (1.0+zdqs(i, k))(v_cptt(i,k)+v_cptj(k)) + RLVTT_generic(local_q(i,k)-zqs(i,k)) )
	239	v_ssig = v_ssig + zdp *(1.0+zdqs(i,k))
	240	ENDDO
	241
	242
	243	! this right here is where the adjustment is done???
	244	DO k = k1, k2
	245	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
	246	cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
	247	v_cptt(i,k)=cp_new_t(k)
	248	local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT_generic
	249	local_t(i,k) = cp_new_t(k) / RCPD
	250
	251	v_t = MAX(local_t(i,k),15.)
	252	v_p = pplay(i,k)
	253
	254	call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k))
	255	call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
	256
	257	!call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
	258	!call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
	259
	260	ENDDO
	261	Enddo ! nn=1,niter
	262
	263
	264	!--------------------------------------------------- sounding downwards
	265	! -- we refine the prognostic variables in
	266	! -- the layer about to be adjusted
	267
	268	! DO k = k1, k2
	269	! v_cptt(i,k) = RCPD * local_t(i,k)
	270	! v_t = local_t(i,k)
	271	! v_p = pplay(i,k)
	272	!
	273	! call Psat_water(v_t,v_p,zpsat,zqs(i,k))
	274	! call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
	275	! ENDDO
	276
	277	DO k = 2, nlayer
	278	zdpm = pplev(i,k-1) - pplev(i,k)
	279	zdp = pplev(i,k) - pplev(i,k+1)
	280	! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
	281	! + RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
	282	! * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	283	! / (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
	284	! general case where water is not a trace gas (JL13)
	285	v_zqs = (zqs(i,k-1)zdpm + zqs(i,k)zdp)/(zdpm+zdp)
	286	v_cptt2 = (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp)/(zdpm+zdp)
	287	v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))zdpm + (1.-zpsat(i,k)/pplay(i,k))zdp)/(zdpm+zdp)
	288	v_dlnpsat = (zdlnpsat(i,k-1)zdpm + zdlnpsat(i,k)zdp)/(zdpm+zdp)
	289	gamcpdz(i,k) = ( (R/RCPDv_cptt2(1.- v_zqs) + RLVTT_genericv_zqs) (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) &
	290	/ (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)v_pratio + v_zqs v_dlnpsat)
	291	ENDDO
	292
	293	! Test to see if we've reached the bottom
	294
	295	IF (k1 .EQ. 1) GOTO 841 ! yes we have!
	296	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
	297	zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) &
	298	+ (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
	299	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
	300
	301	840 CONTINUE
	302	k1 = k1 - 1
	303	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	304	zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1)) &
	305	*(pplev(i,k1-1)-pplev(i,k1))
	306	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
	307	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
	308	GOTO 840
	309	ELSE
	310	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	311	ENDIF
	312	841 CONTINUE
	313
	314	GOTO 810 ! look for other layers higher up
	315
	316	9999 CONTINUE ! loop over all the points
	317
	318	!-----------------------------------------------------------------------
	319	! Determine the cloud fraction (hypothese: la nebulosite a lieu
	320	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
	321
	322	DO k = 1, nlayer
	323	DO i = 1, ngrid
	324	IF (itest(i)) THEN
	325	delta_q(i,k) = local_q(i,k) - zq(i,k)
	326	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
	327	ENDIF
	328	ENDDO
	329	ENDDO
	330
	331	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
	332	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
	333	! diminue et d'une maniere proportionnelle a cet diminution):
	334	DO i = 1, ngrid
	335	IF (itest(i)) THEN
	336	zq1(i) = 0.0
	337	zq2(i) = 0.0
	338	ENDIF
	339	ENDDO
	340	DO k = 1, nlayer
	341	DO i = 1, ngrid
	342	IF (itest(i)) THEN
	343	zdp = pplev(i,k)-pplev(i,k+1)
	344	zq1(i) = zq1(i) - delta_q(i,k) * zdp
	345	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
	346	ENDIF
	347	ENDDO
	348	ENDDO
	349	DO k = 1, nlayer
	350	DO i = 1, ngrid
	351	IF (itest(i)) THEN
	352	IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
	353	ENDIF
	354	ENDDO
	355	ENDDO
	356
	357	DO k = 1, nlayer
	358	DO i = 1, ngrid
	359	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
	360	ENDDO
	361	ENDDO
	362
	363	DO k = 1, nlayer
	364	DO i = 1, ngrid
	365	d_t(i,k) = local_t(i,k) - zt(i,k)
	366	d_q(i,k) = local_q(i,k) - zq(i,k)
	367	ENDDO
	368	ENDDO
	369
	370	! now subroutine -----> GCM variables
	371	DO k = 1, nlayer
	372	DO i = 1, ngrid
	373
	374	pdtmana(i,k) = d_t(i,k)/ptimestep
	375	pdqmana(i,k,i_vap_generic) = d_q(i,k)/ptimestep
	376	pdqmana(i,k,i_ice_generic) = d_ql(i,k)/ptimestep
	377
	378	ENDDO
	379	ENDDO
	380
	381	end subroutine moistadj_generic

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

Download in other formats:

Original Format