Context Navigation

condensation_generic_mod.F90 @ 3058

Last change on this file since 3058 was 3044, checked in by jleconte, 15 months ago
perfect_vap_profile option
File size: 14.9 KB

Rev	Line
[2701]	1	module condensation_generic_mod
[2690]	2	implicit none
	3
	4	contains
	5
[2701]	6	subroutine condensation_generic(ngrid,nlayer,nq,ptimestep, pplev, pplay, &
[2724]	7	pt, pq, pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
[2690]	8	use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity
[2711]	9	use generic_cloud_common_h
[2690]	10	USE tracer_h
[2717]	11	USE mod_phys_lmdz_para, only: is_master
[2720]	12	use generic_tracer_index_mod, only: generic_tracer_index
[2690]	13	IMPLICIT none
	14
	15	!=======================================================================
	16	!
	17	! Purpose
	18	! -------
	19	! Calculates large-scale condensation of generic tracer "tname".
	20	! By convention, tname ends with the suffix "_vap", as it represents the
	21	! gas phase of the generic tracer
	22	!
	23	! Authors
	24	! -------
[2714]	25	! Adapted from largescale.F90 by Lucas Teinturier & Noé Clément (2022)
[2690]	26	! largescale.F90 adapted from the LMDTERRE code by R. Wordsworth (2009)
	27	! Original author Z. X. Li (1993)
	28	!
	29	!=========================================================================
	30
	31	INTEGER, intent(in) :: ngrid,nlayer,nq
	32
	33	! Arguments
	34	REAL, intent(in) :: ptimestep ! intervalle du temps (s)
	35	REAL, intent(in) :: pplev(ngrid,nlayer+1) ! pression a inter-couche
	36	REAL, intent(in) :: pplay(ngrid,nlayer) ! pression au milieu de couche
	37	REAL, intent(in) :: pt(ngrid,nlayer) ! temperature (K)
	38	REAL, intent(in) :: pq(ngrid,nlayer,nq) ! tracer mixing ratio (kg/kg)
	39	REAL, intent(in) :: pdt(ngrid,nlayer) ! physical temperature tendency (K/s)
	40	REAL, intent(in) :: pdq(ngrid,nlayer,nq) ! physical tracer tendency (K/s)
[2700]	41	! CHARACTER(*), intent(in) :: tname_vap ! name of the tracer we consider. BY convention, it ends with _vap !!!
[2690]	42	REAL, intent(out) :: pdtlsc(ngrid,nlayer) ! incrementation de la temperature (K)
[2700]	43	REAL, intent(out) :: pdqvaplsc(ngrid,nlayer,nq) ! incrementation de la vapeur du traceur
	44	REAL, intent(out) :: pdqliqlsc(ngrid,nlayer,nq) ! incrementation du traceur liquide
[2724]	45	REAL, intent(out) :: rneb(ngrid,nlayer,nq) ! fraction nuageuse
[2690]	46
	47	! Options :
	48	real, save :: metallicity !metallicity of planet
[3044]	49	REAL, SAVE :: qvap_deep ! deep mixing ratio of vapor when simulating bottom less planets
	50	REAL, SAVE :: qvap_top ! top mixing ratio of vapor when simulating bottom less planets
	51	logical, save :: perfect_vap_profile
	52	!$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, perfect_vap_profile)
[2690]	53
	54	! Local variables
[2720]	55
	56	! to call only one time the ice/vap pair of a tracer
[2722]	57	logical call_ice_vap_generic
[2720]	58
[2690]	59	INTEGER i, k , nn, iq
	60	INTEGER,PARAMETER :: nitermax=5000
[3044]	61	REAL tau ! tau is in seconds and must not be lower than the physical step time.
	62	integer k_cold_trap
[2690]	63	DOUBLE PRECISION,PARAMETER :: alpha=.1,qthreshold=1.d-8
	64	! JL13: if "careful, T<Tmin in psat water" appears often, you may want to stabilise the model by
	65	! decreasing alpha and increasing nitermax accordingly
	66	DOUBLE PRECISION zq(ngrid)
	67	DOUBLE PRECISION zcond(ngrid),zcond_iter
	68	DOUBLE PRECISION zqs(ngrid)
	69	real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,Lcp,zqs_temp,zdqs
[3044]	70	real zqs_temp_1, zqs_temp_2, zqs_temp_3
[2690]	71	integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
[2700]	72	! CHARACTER(len=*) :: tname_ice
[2690]	73	! evaporation calculations
	74	REAL dqevap(ngrid,nlayer),dtevap(ngrid,nlayer)
	75	REAL qevap(ngrid,nlayer,nq)
	76	REAL tevap(ngrid,nlayer)
	77
	78	DOUBLE PRECISION zx_q(ngrid)
[2890]	79	DOUBLE PRECISION zy_q(ngrid)
[2690]	80	LOGICAL,SAVE :: firstcall=.true.
	81	!$OMP THREADPRIVATE(firstcall)
	82	IF (firstcall) THEN
	83	write(,) "value for metallicity? "
	84	metallicity=0.0 ! default value
	85	call getin_p("metallicity",metallicity)
	86	write(,) " metallicity = ",metallicity
[2720]	87
[3044]	88	write(,) "Deep generic tracer vapor mixing ratio ? (no effect if negative) "
[2720]	89	qvap_deep=-1. ! default value
	90	call getin_p("qvap_deep",qvap_deep)
	91	write(,) " qvap_deep = ",qvap_deep
	92
[3044]	93	write(,) "top generic tracer vapor mixing ratio ? (no effect if negative) "
	94	qvap_top=-1. ! default value
	95	call getin_p("qvap_top",qvap_top)
	96	write(,) " qvap_top = ",qvap_top
	97
	98	write(,) " perfect_vap_profile ? "
	99	perfect_vap_profile=.false. ! default value
	100	call getin_p("perfect_vap_profile",perfect_vap_profile)
	101	write(,) " perfect_vap_profile = ",perfect_vap_profile
	102
[2690]	103	firstcall = .false.
	104	ENDIF
[2700]	105	! Initialisation of outputs and local variables
[2690]	106	pdtlsc(1:ngrid,1:nlayer) = 0.0
[2700]	107	pdqvaplsc(1:ngrid,1:nlayer,1:nq) = 0.0
	108	pdqliqlsc(1:ngrid,1:nlayer,1:nq) = 0.0
	109	dqevap(1:ngrid,1:nlayer)=0.0
	110	dtevap(1:ngrid,1:nlayer)=0.0
	111	qevap(1:ngrid,1:nlayer,1:nq)=0.0
	112	tevap(1:ngrid,1:nlayer)=0.0
[2724]	113	rneb(1:ngrid,1:nlayer,1:nq) = 0.0
[2700]	114	! Let's loop on tracers
[2690]	115	do iq=1,nq
[2706]	116
[2722]	117	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
[2706]	118
[2722]	119	if(call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
[2711]	120	m=constants_mass(iq)
	121	delta_vapH=constants_delta_vapH(iq)
	122	Tref=constants_Tref(iq)
	123	Pref=constants_Pref(iq)
[2725]	124	epsi_generic=constants_epsi_generic(iq)
	125	RLVTT_generic=constants_RLVTT_generic(iq)
[2711]	126	metallicity_coeff=constants_metallicity_coeff(iq)
	127
[2725]	128	Lcp=RLVTT_generic/cpp ! need to be init here
[2690]	129
[2700]	130	! Vertical loop (from top to bottom)
	131	DO k = nlayer, 1, -1
[2704]	132	zt(1:ngrid)=pt(1:ngrid,k)+pdt(1:ngrid,k)*ptimestep
[2690]	133
[2700]	134	! Computes Psat and the partial condensation
	135	DO i = 1, ngrid
	136	local_p=pplay(i,k)
	137	if(zt(i).le.15.) then
	138	print*,'in lsc',i,k,zt(i)
	139	! zt(i)=15. ! check too low temperatures
	140	endif
	141	zx_q(i) = pq(i,k,igcm_generic_vap)+pdq(i,k,igcm_generic_vap)*ptimestep
[2724]	142
[2890]	143	call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp)
	144	zy_q(i) = pq(i,k,igcm_generic_ice)+pdq(i,k,igcm_generic_ice)*ptimestep
	145
	146	if ((zx_q(i) .le. zqs_temp) .and. (zy_q(i) .eq. 0.)) then
	147	! if we are are not saturated and if there is no ice
	148	! then no change
	149
	150	zcond(i) = 0.0d0
	151
	152	else ! if we are saturated : we must evaporate
	153	! if there is ice : we must check if we can condensate
	154
	155	! iterative process to stabilize the scheme when large water amounts JL12
	156	zcond(i) = 0.0d0
	157	Do nn=1,nitermax
	158	call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp)
	159	zqs(i)=zqs_temp
	160	call Lcpdqsat_generic(zt(i),local_p,psat_tmp,zqs_temp,zdqs,dlnpsat_tmp)
	161	zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs)
	162	!zcond can be negative here
	163	zx_q(i) = zx_q(i) - zcond_iter
	164	zcond(i) = zcond(i) + zcond_iter
	165	zt(i) = zt(i) + zcond_iter*Lcp
	166	if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit
	167	if (nn.eq.nitermax) print*,'itermax in largescale'
	168	End do ! niter
	169
	170	! if zcond(i) > 0, zcond(i) is the amount of vapor that we can condensate
	171	! because we are saturated. zcond(i) will not change below
	172	! if zcond(i) < 0, zcond(i) is the amount of ice that we can evaporate.
	173	! We can not evaporate more than the existing ice,
	174	! so the line below is to check how much we can evaporate.
	175	! If there is no ice available, zcond(i) will be 0. (first condidition of "if")
	176
	177	zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_generic_ice)+pdq(i,k,igcm_generic_ice)*ptimestep))
	178
	179	endif
	180
[2724]	181	if (zcond(i) .gt. 0.) then
	182	rneb(i,k,iq)=1
	183	else
	184	rneb(i,k,iq)=0.
	185	endif
	186
[2700]	187	zcond(i) = zcond(i)/ptimestep
	188	ENDDO ! i=1,ngrid
[2690]	189
[2700]	190	!Tendances de t et q
[2704]	191	pdqvaplsc(1:ngrid,k,igcm_generic_vap) = - zcond(1:ngrid)
[2700]	192	pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_vap)
	193	pdtlsc(1:ngrid,k) = pdtlsc(1:ngrid,k) + pdqliqlsc(1:ngrid,k,igcm_generic_ice)*Lcp
[2690]	194
[2700]	195	Enddo ! k= nlayer, 1, -1
[2720]	196
[3044]	197	if ((perfect_vap_profile) .and. (ngrid.eq.1)) then
	198	! perfect_vap_profile is a mode that should a priori only be used in 1D:
	199	! it aims to force the vap profile:
	200	! - below condensation, it is forced to qvap_deep
	201	! - at condensation levels, it is forced to 99% of sat
	202	! - above the cold trap, it is forced to the value allowed by the cold trap
	203	tau=3243600
	204
	205	k_cold_trap = 2
	206	do k=2,nlayer-1
	207
	208	zt(1)=pt(1,k-1)+pdt(1,k-1)*ptimestep
	209	call Psat_generic(zt(1),pplay(1,k-1),metallicity,psat_tmp,zqs_temp_1)
	210	zt(1)=pt(1,k)+pdt(1,k)*ptimestep
	211	call Psat_generic(zt(1),pplay(1,k),metallicity,psat_tmp,zqs_temp_2)
	212	zt(1)=pt(1,k+1)+pdt(1,k+1)*ptimestep
	213	call Psat_generic(zt(1),pplay(1,k+1),metallicity,psat_tmp,zqs_temp_3)
	214
	215	if ((zqs_temp_1 .gt. zqs_temp_2) .and. (zqs_temp_3 .gt. zqs_temp_2)) then
	216	k_cold_trap = k
	217	exit
	218	endif
	219	enddo
	220	if (k_cold_trap .lt. nlayer) then
	221	do k=k_cold_trap+1,nlayer
	222	pdqvaplsc(1,k,igcm_generic_vap) = (pq(1,k_cold_trap,igcm_generic_vap) - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)
	223	enddo
	224	endif
	225
	226	do k=1,k_cold_trap
	227	zt(1)=pt(1,k)+pdt(1,k)*ptimestep
	228	call Psat_generic(zt(1),pplay(1,k),metallicity,psat_tmp,zqs_temp)
	229	if (zqs_temp .gt. qvap_deep) then
	230	pdqvaplsc(1,k,igcm_generic_vap) = (qvap_deep - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)
	231	endif
	232	if (zqs_temp .lt. qvap_deep) then
	233	pdqvaplsc(1,k,igcm_generic_vap) = (0.99*zqs_temp - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)
	234	endif
	235	enddo
	236
	237	pdqliqlsc(1:ngrid,:,igcm_generic_ice) = 0.
	238	pdtlsc(1:ngrid,:) = 0.
	239	endif
	240
[2720]	241	if (qvap_deep >= 0.) then
[3044]	242	if (ngrid.eq.1) then ! if 1D
	243	tau=3243600 ! tau must not be lower than the physical step time. In 1D time step is very long
	244	else
	245	tau=3600 ! for 3D
	246	endif
[2720]	247	!brings lower generic vapor ratio to a fixed value.
[2890]	248	! tau is in seconds and must not be lower than the physical step time.
[2720]	249	pdqvaplsc(1:ngrid,1,igcm_generic_vap) = (qvap_deep - pq(1:ngrid,1,igcm_generic_vap))/tau - pdq(1:ngrid,1,igcm_generic_vap)
	250	endif
[3044]	251	if (qvap_top >= 0.) then
	252	if (ngrid.eq.1) then ! if 1D
	253	tau=3243600 ! tau must not be lower than the physical step time. In 1D time step is very long
	254	else
	255	tau=3600 ! for 3D
	256	endif
	257	!brings lower generic vapor ratio to a fixed value.
	258	! tau is in seconds and must not be lower than the physical step time.
	259	pdqvaplsc(1:ngrid,nlayer,igcm_generic_vap) = (qvap_top - pq(1:ngrid,nlayer,igcm_generic_vap))/tau - pdq(1:ngrid,nlayer,igcm_generic_vap)
	260	endif
[2722]	261	endif !if(call_ice_vap_generic)
[2700]	262	enddo ! iq=1,nq
[2701]	263
	264	end subroutine condensation_generic
	265	end module condensation_generic_mod

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

Download in other formats:

Original Format