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

Line
1	module condensation_generic_mod
2	implicit none
3
4	contains
5
6	subroutine condensation_generic(ngrid,nlayer,nq,ptimestep, pplev, pplay, &
7	pt, pq, pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
8	use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity
9	use generic_cloud_common_h
10	USE tracer_h
11	USE mod_phys_lmdz_para, only: is_master
12	use generic_tracer_index_mod, only: generic_tracer_index
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	! -------
25	! Adapted from largescale.F90 by Lucas Teinturier & Noé Clément (2022)
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)
41	! CHARACTER(*), intent(in) :: tname_vap ! name of the tracer we consider. BY convention, it ends with _vap !!!
42	REAL, intent(out) :: pdtlsc(ngrid,nlayer) ! incrementation de la temperature (K)
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
45	REAL, intent(out) :: rneb(ngrid,nlayer,nq) ! fraction nuageuse
46
47	! Options :
48	real, save :: metallicity !metallicity of planet
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)
53
54	! Local variables
55
56	! to call only one time the ice/vap pair of a tracer
57	logical call_ice_vap_generic
58
59	INTEGER i, k , nn, iq
60	INTEGER,PARAMETER :: nitermax=5000
61	REAL tau ! tau is in seconds and must not be lower than the physical step time.
62	integer k_cold_trap
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
70	real zqs_temp_1, zqs_temp_2, zqs_temp_3
71	integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
72	! CHARACTER(len=*) :: tname_ice
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)
79	DOUBLE PRECISION zy_q(ngrid)
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
87
88	write(,) "Deep generic tracer vapor mixing ratio ? (no effect if negative) "
89	qvap_deep=-1. ! default value
90	call getin_p("qvap_deep",qvap_deep)
91	write(,) " qvap_deep = ",qvap_deep
92
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
103	firstcall = .false.
104	ENDIF
105	! Initialisation of outputs and local variables
106	pdtlsc(1:ngrid,1:nlayer) = 0.0
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
113	rneb(1:ngrid,1:nlayer,1:nq) = 0.0
114	! Let's loop on tracers
115	do iq=1,nq
116
117	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
118
119	if(call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
120	m=constants_mass(iq)
121	delta_vapH=constants_delta_vapH(iq)
122	Tref=constants_Tref(iq)
123	Pref=constants_Pref(iq)
124	epsi_generic=constants_epsi_generic(iq)
125	RLVTT_generic=constants_RLVTT_generic(iq)
126	metallicity_coeff=constants_metallicity_coeff(iq)
127
128	Lcp=RLVTT_generic/cpp ! need to be init here
129
130	! Vertical loop (from top to bottom)
131	DO k = nlayer, 1, -1
132	zt(1:ngrid)=pt(1:ngrid,k)+pdt(1:ngrid,k)*ptimestep
133
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
142
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
181	if (zcond(i) .gt. 0.) then
182	rneb(i,k,iq)=1
183	else
184	rneb(i,k,iq)=0.
185	endif
186
187	zcond(i) = zcond(i)/ptimestep
188	ENDDO ! i=1,ngrid
189
190	!Tendances de t et q
191	pdqvaplsc(1:ngrid,k,igcm_generic_vap) = - zcond(1:ngrid)
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
194
195	Enddo ! k= nlayer, 1, -1
196
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
241	if (qvap_deep >= 0.) then
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
247	!brings lower generic vapor ratio to a fixed value.
248	! tau is in seconds and must not be lower than the physical step time.
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
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
261	endif !if(call_ice_vap_generic)
262	enddo ! iq=1,nq
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