Context Navigation

moistadj_generic.F90 @ 3436

Last change on this file since 3436 was 3279, checked in by jleconte, 8 months ago
flag moist_convection_inhibition (Noe C)
File size: 15.4 KB

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

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

Download in other formats:

Original Format