Context Navigation

co2condens_mod4micro.F90 @ 2566

Last change on this file since 2566 was 2524, checked in by cmathe, 4 years ago
minor fix in co2condens4micro; fix co2 sedimentation rate outputs
File size: 35.5 KB

Line
1	!======================================================================================================================!
2	! Module: CO2 condensation for the CO2 cloud microphysics =============================================================!
3	!----------------------------------------------------------------------------------------------------------------------!
4	! Authors: Christophe Mathé, Anni Määttänen
5	! Date: 16/04/2020
6	!----------------------------------------------------------------------------------------------------------------------!
7	! Contains subroutines:
8	! - co2condens4micro: condensation/sublimation of CO2 ice on the ground and compute pressure change resulting
9	!
10	! - vl1d: Van-Leer scheme
11	!======================================================================================================================!
12	module co2condens_mod4micro
13
14	implicit none
15
16	contains
17	!======================================================================================================================!
18	! SUBROUTINE: co2condens4micro ========================================================================================!
19	!======================================================================================================================!
20	! Subject:
21	!---------
22	! Condensation/sublimation of CO2 ice on the ground and compute pressure change resulting
23	!----------------------------------------------------------------------------------------------------------------------!
24	! Comments:
25	!----------
26	! Adapted from co2condens_mod.F
27	!----------------------------------------------------------------------------------------------------------------------!
28	! Paper:
29	!-------
30	! Forget et al. (2008), "Non condensable gas enrichment and depletion in the Martian polar regions."
31	!----------------------------------------------------------------------------------------------------------------------!
32	! Algorithm:
33	!-----------
34	! 1. Initialization
35	! 2. Firstcall
36	! 3. Compute CO2 Volume mixing ratio
37	! 4. Set zcondicea, zfallice from co2clouds condensation rate and set zt
38	! 5. Main co2condens
39	! 5.1. Forecast of ground temperature ztsrf and frost temperature ztcondsol
40	! 5.2. Check if we have condensation/sublimation on the ground
41	! 5.3. Compute zfallheat
42	! 5.4. Compute direct condensation/sublimation of CO2 ice
43	! 5.4.a. If there is not enough CO2 tracer in 1st layer to condense
44	! 5.4.b. If the entire CO2 ice layer sublimes (including what has just condensed in the atmosphere)
45	! 5.5. Changing CO2 ice amount and pressure
46	! 5.6. Surface albedo and emissivity of the surface below the snow (emisref)
47	! 5.6.a. Check that amont of CO2 ice is not problematic
48	! 5.6.b. Set albedo and emissivity of the surface
49	! 5.6.c. Set pemisurf to emissiv when there is bare surface (needed for co2snow)
50	! 5.7. Correction to account for redistribution between sigma or hybrid layers when changing surface pressure (and
51	! warming/cooling of the CO2 currently changing phase).
52	! 5.7.a. Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
53	! 5.7.b. Mass of each layer at the end of timestep
54	! 5.7.c. Corresponding fluxes for T, U, V and Q (averaging operator for transport)
55	! 5.7.c.i. Value transfert at the surface interface when condensation/sublimation
56	! 5.7.c.ii. Van Leer scheme
57	! 5.7.c.iii Compute tendencies on T, U, V, Q
58	! 6. CO2 snow / clouds scheme
59	! 7. Extra special case for surface temperature tendency pdtsrfc
60	!----------------------------------------------------------------------------------------------------------------------!
61	subroutine co2condens4micro(ngrid, nlayer, nq, ptimestep, pcapcal, pplay, pplev, ptsrf, pt, pphi, pdt, pdu, pdv, &
62	pdtsrf, pu, pv, pq, pdq, piceco2, psolaralb, pemisurf, pdtc, pdtsrfc, pdpsrf, pduc, &
63	pdvc, pdqc, fluxsurf_sw, zls, zdqssed_co2, pcondicea_co2microp)
64
65	use tracer_mod, only: noms, igcm_co2, igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice, igcm_dust_number, igcm_dust_mass, &
66	igcm_ccnco2_number, igcm_ccnco2_mass, igcm_ccn_mass
67
68	use surfdat_h, only: emissiv, phisfi
69
70	use geometry_mod, only: latitude, longitude_deg, latitude_deg
71
72	use planete_h, only: obliquit
73
74	use comcstfi_h, only: cpp, g, r, pi
75
76	#ifndef MESOSCALE
77	use vertical_layers_mod, only: ap, bp
78	#endif
79
80	implicit none
81
82	include "callkeys.h"
83	!----------------------------------------------------------------------------------------------------------------------!
84	! VARIABLE DECLARATION
85	!----------------------------------------------------------------------------------------------------------------------!
86	! Input arguments:
87	!-----------------
88	integer, intent(in) :: &
89	nq, &! number of tracers
90	ngrid, &! number of atmospheric columns
91	nlayer ! number of vertical layers
92
93	real, intent(in) :: &
94	ptimestep, &! physics timestep (s)
95	pcapcal(ngrid), &! surface specific heat
96	pplay(ngrid,nlayer), &! mid-layer pressure (Pa)
97	pplev(ngrid,nlayer+1), &! inter-layer pressure (Pa)
98	ptsrf(ngrid), &! surface temperature (K)
99	pphi(ngrid,nlayer), &! geopotential (m2.s-2)
100	pt(ngrid,nlayer), &! atmospheric temperature (K)
101	pu(ngrid,nlayer), &! zonal wind (m/s)
102	pv(ngrid,nlayer), &! meridional wind (m/s)
103	pdt(ngrid,nlayer), &! tendency on temperature from previous physical processes (K/s)
104	pdu(ngrid,nlayer), &! tendency on zonal wind from previous physical processes(m/s2)
105	pdv(ngrid,nlayer), &! tendency on meridional wind from previous physical processes (m/s2)
106	pdtsrf(ngrid), &! tendency on surface temperature from previous physical processes (K/s)
107	pq(ngrid,nlayer,nq), &! tracers (../kg_air)
108	pdq(ngrid,nlayer,nq), &! tendency on tracers from previous physical processes
109	zdqssed_co2(ngrid), &! CO2 flux at the surface (kg.m-2.s-1)
110	fluxsurf_sw(ngrid,2), &! added to calculate flux dependent albedo
111	zls, &! solar longitude (rad)
112	pcondicea_co2microp(ngrid,nlayer) ! tendency due to CO2 condensation (kg/kg.s-1)
113	!----------------------------------------------------------------------------------------------------------------------!
114	! In/output arguments:
115	!---------------------
116	real, intent(inout) :: &
117	piceco2(ngrid), &! CO2 ice on the surface (kg.m-2)
118	pemisurf(ngrid), &! emissivity of the surface
119	psolaralb(ngrid,2) ! albedo of the surface
120	!----------------------------------------------------------------------------------------------------------------------!
121	! Output arguments:
122	!------------------
123	real, intent(out) :: &
124	pdtc(ngrid,nlayer), &! tendency on temperature dT/dt due to cond/sub (K/s)
125	pdtsrfc(ngrid), &! tendency on surface temperature (K/s)
126	pdpsrf(ngrid), &! tendency on surface pressure (Pa/s)
127	pduc(ngrid,nlayer), &! tendency on zonal wind (m.s-2)
128	pdvc(ngrid,nlayer), &! tendency on meridional wind (m.s-2)
129	pdqc(ngrid,nlayer,nq) ! tendency on tracers
130	!----------------------------------------------------------------------------------------------------------------------!
131	! Local:
132	!-------
133	!----1) Parameters:
134	!------------------
135	real, parameter :: &
136	latcond = 595594, &! latent heat of solid CO2 ice (J/kg)
137	cpice = 1000., &! specific heat of CO2 ice (J.kg-1.K-1)
138	tcond1mb = 136.27,&! condensation temperature at 1 mbar (K)
139	m_co2 = 44.01E-3, &! CO2 molecular mass (kg/mol)
140	m_noco2 = 33.37E-3 ! non condensible molecular mass (kg/mol)
141
142	logical, parameter :: &
143	improved_ztcond = .true. ! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer) then
144	! condensation temperature is computed using partial pressure of CO2
145	!----------------------------------------------------------------------------------------------------------------------!
146	!----2) Saved:
147	!-------------
148	real, save :: &
149	A, &! coefficient used to compute mean molecular mass
150	B, &! coefficient used to compute mean molecular mass
151	acond,&! coefficient used to compute ztcondsol
152	bcond,&! coefficient used to compute ztcondsol
153	ccond ! coefficient used to compute ztcondsol
154
155	logical, save :: &
156	firstcall = .true. ! Used to compute saved variables
157	!----------------------------------------------------------------------------------------------------------------------!
158	!----3) Variables:
159	!-----------------
160	integer :: &
161	l, &! loop on layers
162	ig, &! loop on ngrid points
163	iq ! loop on tracer
164
165	real :: &
166	qco2, &! effective quantity of CO2, used to compute mean molecular mass
167	mmean, &! mean molecular mass
168	zfallheat, &! aerodynamical friction and energy used to heat the amount of ice
169	zmflux(nlayer+1), &! mass fluxes through the sigma levels (kg.m-2.s-1)
170	zu(nlayer), &! effective zonal wind
171	zv(nlayer), &! effective meridional wind
172	zq(nlayer,nq), &! effective tracers quantities
173	zq1(nlayer), &! buffer of zq
174	ztsrf(ngrid), &! effective temperature at the surface
175	ztm(nlayer+1), &! temperature fluxes through the sigma levels
176	zum(nlayer+1), &! zonal wind fluxes through the sigma levels
177	zvm(nlayer+1), &! meridional wind fluxes through the sigma levels
178	zqm(nlayer+1,nq), &! quantity of tracers flux through the sigma levels
179	zqm1(nlayer), &! quantity of tracers after Van-Leer scheme
180	masse(nlayer), &! mass layer (kg)
181	w(nlayer+1), &! total mass fluxes through the sigma levels during ptimestep (kg.m-2)
182	availco2, &! available quantity of co2 (kg)
183	emisref(ngrid), &! emissivity reference
184	vmr_co2(ngrid,nlayer), &! CO2 volume mixing ratio
185	zt(nlayer), &! effective temperature in the atmosphere (K)
186	ztcond(ngrid,nlayer+1), &! CO2 condensation temperature (atm)
187	ztcondsol(ngrid), &! CO2 condensation temperature (surface)
188	zdiceco2(ngrid), &! tendency on co2ice surface tracer (kg/m2/s)
189	zcondicea(ngrid,nlayer),&! condensation rate in layer l (kg/m2/s)
190	zcondices(ngrid), &! condensation rate on the ground (kg/m2/s)
191	zfallice(ngrid) ! amount of ice falling from layer l (kg/m2/s)
192
193	logical :: &
194	condsub(ngrid) ! True if there is condensation/sublimation (used for co2snow)
195
196
197	! check with co2 cloud parameterisation
198	real :: &
199	zt_2(ngrid,nlayer), &
200	ztcond_2(ngrid,nlayer+1), &
201	zfallice_2(ngrid,nlayer+1), &
202	pdtc_2(ngrid,nlayer), &
203	zfallheat_2, &
204	zcondicea_2(ngrid,nlayer), &
205	diff_zcondicea(ngrid,nlayer), &
206	diff_zfallice(ngrid)
207	!======================================================================================================================!
208	! BEGIN ===============================================================================================================!
209	!======================================================================================================================!
210	! 1. Initialization
211	!----------------------------------------------------------------------------------------------------------------------!
212	availco2 = 0.
213	zfallheat = 0.
214	zt(1:nlayer) = 0.
215	ztcond(1:ngrid, 1:nlayer+1) = 0.
216	ztcondsol(1:ngrid) = 0.
217	zmflux(1:nlayer+1) = 0.
218	zu(1:nlayer) = 0.
219	zv(1:nlayer) = 0.
220	zq(1:nlayer, 1:nq) = 0.
221	zq1(1:nlayer) = 0.
222	ztsrf(1:ngrid) = 0.
223	ztm(1:nlayer+1) = 0.
224	zum(1:nlayer+1) = 0.
225	zvm(1:nlayer+1) = 0.
226	zqm(1:nlayer+1, 1:nq) = 0.
227	masse(1:nlayer) = 0.
228	w(1:nlayer+1) = 0.
229	emisref(1:ngrid) = 0.
230	vmr_co2(1:ngrid, 1:nlayer) = 0.
231	zcondices(1:ngrid) = 0.
232	pdtsrfc(1:ngrid) = 0.
233	pdpsrf(1:ngrid) = 0.
234	zdiceco2(1:ngrid) = 0.
235	condsub(1:ngrid) = .false.
236	zcondicea(1:ngrid, 1:nlayer) = 0.
237	zfallice(1:ngrid) = 0.
238	pduc(1:ngrid, 1:nlayer) = 0.
239	pdvc(1:ngrid, 1:nlayer) = 0.
240	pdqc(1:ngrid, 1:nlayer, 1:nq) = 0.
241	pdtc(1:ngrid,1:nlayer) = 0.
242	!----------------------------------------------------------------------------------------------------------------------!
243	! 2. Firstcall
244	!----------------------------------------------------------------------------------------------------------------------!
245	! AS: firstcall OK absolute
246	if (firstcall) then
247	firstcall = .false.
248
249	bcond = 1. / tcond1mb
250	ccond = cpp / (g*latcond)
251	acond = r / latcond
252
253	write(,)'CO2condens: improved_ztcond=', improved_ztcond
254	write(,)'In co2condens:Tcond(P=1mb)=', tcond1mb, ' Lcond=', latcond
255	write(,)'acond,bcond,ccond', acond, bcond, ccond
256
257	! Prepare Special treatment if one of the tracer is CO2 gas. Compute A and B coefficient use to compute mean molecular
258	! mass Mair defined by:
259	! 1/Mair = q(igcm_co2)/m_co2 + (1-q(igcm_co2))/m_noco2
260	! 1/Mair = A*q(igcm_co2) + B
261	A = (1./m_co2 - 1./m_noco2)
262	B = 1./m_noco2
263	end if ! of IF (firstcall)
264	!----------------------------------------------------------------------------------------------------------------------!
265	! 3. Compute CO2 Volume mixing ratio
266	!----------------------------------------------------------------------------------------------------------------------!
267	if (improved_ztcond.and.(igcm_co2/=0)) then
268	do l = 1, nlayer
269	do ig = 1, ngrid
270	qco2 = pq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)*ptimestep
271	! Mean air molecular mass = 1/(q(igcm_co2)/m_co2 + (1-q(igcm_co2))/m_noco2)
272	mmean = 1. / (A*qco2 +B)
273	vmr_co2(ig,l) = (qco2*mmean) / m_co2
274	end do
275	end do
276	else
277	do l = 1, nlayer
278	do ig = 1, ngrid
279	vmr_co2(ig,l) = 0.95
280	end do
281	end do
282	end if
283	!----------------------------------------------------------------------------------------------------------------------!
284	! 4. Set zcondicea, zfallice from co2clouds condensation rate
285	!----------------------------------------------------------------------------------------------------------------------!
286	do l = nlayer, 1, -1
287	do ig = 1, ngrid
288	zcondicea(ig,l) = pcondicea_co2microp(ig,l) * (pplev(ig,l) - pplev(ig,l+1))/g
289	end do
290	end do
291
292	! Only sedimentation falls on the ground !
293	do ig = 1, ngrid
294	zfallice(ig) = zdqssed_co2(ig)
295	piceco2(ig) = piceco2(ig) + zfallice(ig)*ptimestep
296	end do
297
298	! Compute without microphysics
299	diff_zcondicea(1:ngrid, 1:nlayer) = 0.
300	diff_zfallice(1:ngrid) = 0.
301	do l =1, nlayer
302	do ig = 1, ngrid
303	zt_2(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
304	end do
305	end do
306
307	do l = 1, nlayer
308	do ig = 1, ngrid
309	if (pplay(ig,l).ge.1e-4) then
310	ztcond_2(ig,l) = 1. / (bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
311	else
312	ztcond_2(ig,l) = 0.0 !mars Monica
313	endif
314	end do
315	end do
316
317	ztcond_2(:,nlayer+1)=ztcond_2(:,nlayer)
318	zfallice_2(:,:) = 0.
319	zcondicea_2(:,:) = 0.
320	do l = nlayer , 1, -1
321	do ig = 1, ngrid
322	pdtc_2(ig,l) = 0.
323	if ((zt_2(ig,l).lt.ztcond_2(ig,l)).or.(zfallice_2(ig,l+1).gt.0)) then
324	if (zfallice_2(ig,l+1).gt.0) then
325	zfallheat_2 = zfallice_2(ig,l+1)(pphi(ig,l+1)-pphi(ig,l) + cpice(ztcond_2(ig,l+1)-ztcond_2(ig,l)))/latcond
326	else
327	zfallheat_2 = 0.
328	end if
329	pdtc_2(ig,l) = (ztcond_2(ig,l) - zt_2(ig,l))/ptimestep
330	zcondicea_2(ig,l) = (pplev(ig,l)-pplev(ig,l+1))ccondpdtc_2(ig,l)- zfallheat_2
331	! Case when the ice from above sublimes entirely
332	! """""""""""""""""""""""""""""""""""""""""""""""
333	if (zfallice_2(ig,l+1).lt.- zcondicea_2(ig,l)) then
334	zcondicea_2(ig,l)= -zfallice_2(ig,l+1)
335	end if
336	zfallice_2(ig,l) = zcondicea_2(ig,l)+zfallice_2(ig,l+1)
337	end if
338	end do
339	end do
340	diff_zcondicea(:,:) = zcondicea_2(:,:) - zcondicea(:,:)
341	diff_zfallice(:) = zfallice_2(:,1) - zfallice(:)
342	call writediagfi(ngrid, "diff_zfallice", "sans - avec microphysique", "", 2, diff_zfallice)
343	call writediagfi(ngrid, "diff_zcondicea", "sans - avec microphysique", "", 3, diff_zcondicea)
344	!----------------------------------------------------------------------------------------------------------------------!
345	! 5. Main co2condens
346	!----------------------------------------------------------------------------------------------------------------------!
347	do ig = 1, ngrid
348	!----------------------------------------------------------------------------------------------------------------------!
349	! 5.1. Forecast of ground temperature ztsrf and frost temperature ztcondsol
350	!----------------------------------------------------------------------------------------------------------------------!
351	ztcondsol(ig) = 1. / (bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
352	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
353	!----------------------------------------------------------------------------------------------------------------------!
354	! 5.2. Check if we have condensation/sublimation on the ground
355	!----------------------------------------------------------------------------------------------------------------------!
356	! ground condensation \|\| falling snow \|\| ground sublimation
357	!----------------------------------------------------------------------------------------------------------------------!
358	if ((ztsrf(ig)<ztcondsol(ig)) .OR. (zfallice(ig)/=0.) .OR. ((ztsrf(ig)>ztcondsol(ig)) .AND. (piceco2(ig)/=0.))) then
359	condsub(ig) = .true.
360	!----------------------------------------------------------------------------------------------------------------------!
361	! 5.3. Compute zfallheat
362	!----------------------------------------------------------------------------------------------------------------------!
363	zfallheat = 0.
364	!----------------------------------------------------------------------------------------------------------------------!
365	! 5.4. Compute direct condensation/sublimation of CO2 ice
366	!----------------------------------------------------------------------------------------------------------------------!
367	zcondices(ig) = pcapcal(ig) * (ztcondsol(ig)-ztsrf(ig)) / (latcond*ptimestep) - zfallheat
368	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig)) / ptimestep
369	!----------------------------------------------------------------------------------------------------------------------!
370	! 5.4.a. If there is not enough CO2 tracer in 1st layer to condense
371	!----------------------------------------------------------------------------------------------------------------------!
372	! Available CO2 tracer in layer 1 at end of timestep (kg/m2)
373	#ifndef MESOSCALE
374	availco2 = pq(ig,1,igcm_co2) * ( (ap(1)-ap(2)) + (bp(1)-bp(2)) * (pplev(ig,1)/g - zcondices(ig)*ptimestep) )
375	#else
376	availco2 = pq(ig,1,igcm_co2)
377	PRINT*, "MESOSCALE: CO2 tracer AT FIRST LEVEL IS NOT CORRECTED FROM SIGMA LEVELS"
378	#endif
379	if ( zcondices(ig) * ptimestep>availco2 ) then
380	zcondices(ig) = availco2/ptimestep
381	zdiceco2(ig) = zcondices(ig) + zfallice(ig)
382	pdtsrfc(ig) = (latcond/pcapcal(ig)) * (zcondices(ig)+zfallheat)
383	end if
384	!----------------------------------------------------------------------------------------------------------------------!
385	! 5.4.b. If the entire CO2 ice layer sublimes (including what has just condensed in the atmosphere)
386	!----------------------------------------------------------------------------------------------------------------------!
387	if ( (piceco2(ig)/ptimestep) <= -zcondices(ig) ) then
388	zcondices(ig) = -piceco2(ig)/ptimestep
389	pdtsrfc(ig) = (latcond/pcapcal(ig)) * (zcondices(ig)+zfallheat)
390	end if
391	!----------------------------------------------------------------------------------------------------------------------!
392	! 5.5. Changing CO2 ice amount and pressure
393	!----------------------------------------------------------------------------------------------------------------------!
394	zdiceco2(ig) = zcondices(ig) + zfallice(ig) + sum(zcondicea(ig,:))
395	piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep
396	pdpsrf(ig) = -zdiceco2(ig) * g
397
398	if (abs(pdpsrf(ig)*ptimestep)>pplev(ig,1)) then
399	print *, 'STOP in condens'
400	print *, 'condensing more than total mass'
401	print *, 'Grid point ', ig
402	print *, 'Longitude(degrees): ', longitude_deg(ig)
403	print *, 'Latitude(degrees): ', latitude_deg(ig)
404	print *, 'Ps = ', pplev(ig,1)
405	print *, 'd Ps = ', pdpsrf(ig)
406	call abort_physic('co2condens4micro', 'condensing more than total mass', 1)
407	end if
408	!----------------------------------------------------------------------------------------------------------------------!
409	! 5.6. Surface albedo and emissivity of the surface below the snow (emisref)
410	!----------------------------------------------------------------------------------------------------------------------!
411	! 5.6.a. Check that amont of CO2 ice is not problematic
412	!----------------------------------------------------------------------------------------------------------------------!
413	if(.not.piceco2(ig)>=0.) then
414	if(piceco2(ig)<=-5.e-8) then
415	write(,)'WARNING co2condens piceco2(', ig, ')=', piceco2(ig)
416	piceco2(ig) = 0.
417	end if
418	end if
419	!----------------------------------------------------------------------------------------------------------------------!
420	! 5.6.c. Set pemisurf to emissiv when there is bare surface (needed for co2snow)
421	!----------------------------------------------------------------------------------------------------------------------!
422	if (piceco2(ig)==0) then
423	pemisurf(ig) = emissiv
424	end if
425	!----------------------------------------------------------------------------------------------------------------------!
426	! 5.7. Correction to account for redistribution between sigma or hybrid layers when changing surface pressure (and
427	! warming/cooling of the CO2 currently changing phase).
428	!----------------------------------------------------------------------------------------------------------------------!
429	do l= 1, nlayer
430	zt(l) = pt(ig,l) + pdt(ig,l)*ptimestep
431	zu(l) = pu(ig,l) + pdu(ig,l)*ptimestep
432	zv(l) = pv(ig,l) + pdv(ig,l)*ptimestep
433	do iq=1,nq
434	zq(l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
435	end do
436	end do
437	!----------------------------------------------------------------------------------------------------------------------!
438	! 5.7.a. Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
439	!----------------------------------------------------------------------------------------------------------------------!
440	zmflux(1) = - zcondices(ig) - zfallice(ig)
441	do l = 1, nlayer
442	#ifndef MESOSCALE
443	zmflux(l+1) = zmflux(l) - zcondicea(ig,l) &
444	+ (bp(l)-bp(l+1)) * (-pdpsrf(ig)/g)
445	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
446	if (abs(zmflux(l+1))<1E-13.OR.bp(l+1)==0.) then
447	zmflux(l+1) = 0.
448	end if
449	#else
450	zmflux(l+1) = zmflux(l) - zcondicea(ig,l)
451	if (abs(zmflux(l+1))<1E-13) then
452	zmflux(l+1) = 0.
453	end if
454	PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from dPS HAVE TO BE IMPLEMENTED"
455	#endif
456	end do
457	!----------------------------------------------------------------------------------------------------------------------!
458	! 5.7.b. Mass of each layer at the end of timestep
459	!----------------------------------------------------------------------------------------------------------------------!
460	do l = 1, nlayer
461	#ifndef MESOSCALE
462	masse(l) = (pplev(ig,l) - pplev(ig,l+1) + (bp(l)-bp(l+1))pdpsrf(ig)ptimestep)/g
463	#else
464	masse(l) = (pplev(ig,l) - pplev(ig,l+1))/g
465	PRINT*, "MESOSCALE: MASS OF EACH LAYER IS NOT CORRECTED AT END OF TIMESTEP (from SIGMA LEVELS and dPS)"
466	#endif
467	end do
468	!----------------------------------------------------------------------------------------------------------------------!
469	! 5.7.c. Corresponding fluxes for T, U, V and Q (averaging operator for transport)
470	!----------------------------------------------------------------------------------------------------------------------!
471	! 5.7.c.i. Value transfert at the surface interface when condensation/sublimation
472	!----------------------------------------------------------------------------------------------------------------------!
473	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
474	zum(1) = 0
475	zvm(1) = 0
476
477	! Most tracer do not condense
478	do iq = 1, nq
479	zqm(1,iq) = 0.
480	end do
481
482	! Special case if one of the tracer is CO2 gas
483	if (igcm_co2/=0) then
484	zqm(1,igcm_co2) = 1. ! flux is 100% CO2
485	end if
486	!----------------------------------------------------------------------------------------------------------------------!
487	! 5.7.c.ii. Van Leer scheme
488	!----------------------------------------------------------------------------------------------------------------------!
489	do l=1,nlayer
490	w(l)=-zmflux(l)*ptimestep
491	end do
492
493	call vl1d(nlayer,zt,2.,masse,w,ztm)
494	call vl1d(nlayer,zu ,2.,masse,w,zum)
495	call vl1d(nlayer,zv ,2.,masse,w,zvm)
496
497	do iq=1, nq
498	do l=1, nlayer
499	zq1(l) = zq(l,iq)
500	end do
501
502	zqm1(1) = zqm(1,iq)
503	zqm1(2:nlayer) = 0.
504
505	call vl1d(nlayer,zq1,2.,masse,w,zqm1)
506
507	do l = 2, nlayer
508	zq(l,iq) = zq1(l)
509	zqm(l,iq) = zqm1(l)
510	end do
511	end do
512
513	! Surface condensation affects low winds
514	if (zmflux(1)<0) then
515	zum(1) = zu(1) * (w(1)/masse(1))
516	zvm(1) = zv(1) * (w(1)/masse(1))
517	if (w(1)>masse(1)) then ! ensure numerical stability
518	zum(1) = ((zu(1)-zum(2))*masse(1)/w(1)) + zum(2)
519	zvm(1) = ((zv(1)-zvm(2))*masse(1)/w(1)) + zvm(2)
520	end if
521	end if
522
523	ztm(nlayer+1) = zt(nlayer) ! should not be used, but...
524	zum(nlayer+1) = zu(nlayer) ! should not be used, but...
525	zvm(nlayer+1) = zv(nlayer) ! should not be used, but...
526
527	do iq = 1, nq
528	zqm(nlayer+1,iq) = zq(nlayer,iq)
529	end do
530	!----------------------------------------------------------------------------------------------------------------------!
531	! 5.7.c.iii Compute tendencies on T, U, V, Q
532	!----------------------------------------------------------------------------------------------------------------------!
533	#ifdef MESOSCALE
534	! AS: This part must be commented in the mesoscale model
535	! AS: ... to avoid instabilities.
536	! AS: you have to compile with -DMESOSCALE to do so
537	#else
538	do l = 1, nlayer
539	! Tendencies on T
540	pdtc(ig,l) = (1./masse(l)) * ( zmflux(l)(ztm(l) - zt(l)) - zmflux(l+1)(ztm(l+1) - zt(l)) )
541
542	! Tendencies on U
543	pduc(ig,l) = (1./masse(l)) * ( zmflux(l)(zum(l) - zu(l)) - zmflux(l+1)(zum(l+1) - zu(l)) )
544
545	! Tendencies on V
546	pdvc(ig,l) = (1./masse(l)) * ( zmflux(l)(zvm(l) - zv(l)) - zmflux(l+1)(zvm(l+1) - zv(l)) )
547	end do
548	#endif
549	! Tendencies on Q
550	do iq = 1, nq
551	if (iq==igcm_co2) then
552	do l = 1, nlayer
553	pdqc(ig,l,iq) = (1./masse(l)) * (zmflux(l)(zqm(l,iq) - zq(l,iq))- zmflux(l+1)(zqm(l+1,iq) - zq(l,iq))&
554	+ zcondicea(ig,l)*(zq(l,iq) - 1.))
555	end do
556	else
557	do l = 1, nlayer
558	pdqc(ig,l,iq) = (1./masse(l)) * ( zmflux(l)(zqm(l,iq)-zq(l,iq)) - zmflux(l+1)(zqm(l+1,iq)-zq(l,iq))&
559	+ zcondicea(ig,l)*zq(l,iq) )
560	end do
561	end if
562	end do
563	end if ! if
564	end do ! loop on ig
565	!----------------------------------------------------------------------------------------------------------------------!
566	! 5.6.b. Set albedo and emissivity of the surface
567	!----------------------------------------------------------------------------------------------------------------------!
568	call albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
569	!----------------------------------------------------------------------------------------------------------------------!
570	! 6. CO2 snow / clouds scheme
571	!----------------------------------------------------------------------------------------------------------------------!
572	call co2snow(ngrid, nlayer, ptimestep, emisref, condsub, pplev, zcondicea, zcondices, zdqssed_co2, &
573	pemisurf)
574	!----------------------------------------------------------------------------------------------------------------------!
575	! 7. Extra special case for surface temperature tendency pdtsrfc:
576	! We want to fix the south pole temperature to CO2 condensation temperature.
577	!----------------------------------------------------------------------------------------------------------------------!
578	#ifndef MESOSCALE
579	if (caps.and.(obliquit<27.)) then
580	! check if last grid point is the south pole
581	if (abs(latitude(ngrid)-(-pi/2.))<1.e-5) then
582	! NB: Updated surface pressure, at grid point 'ngrid', is ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
583	! write(,)"co2condens: South pole: latitude(ngrid)=", latitude(ngrid)
584	ztcondsol(ngrid) = 1./(bcond-acondlog(.01vmr_co2(ngrid,1) * (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
585	pdtsrfc(ngrid) = (ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
586	end if
587	end if
588	#endif
589	!======================================================================================================================!
590	! END =================================================================================================================!
591	!======================================================================================================================!
592	end subroutine co2condens4micro
593
594
595	!**********************************************************************************************************************!
596	!**********************************************************************************************************************!
597
598
599	!======================================================================================================================!
600	! SUBROUTINE: Van-Leer scheme =========================================================================================!
601	!======================================================================================================================!
602	! Subject:
603	!---------
604	! Operateur de moyenne inter-couche pour calcul de transport type Van-Leer " pseudo amont " dans la verticale
605	!----------------------------------------------------------------------------------------------------------------------!
606	! Comments:
607	!----------
608	! q,w are input arguments for the s-pg ....
609	!----------------------------------------------------------------------------------------------------------------------!
610	! Paper:
611	!-------
612	! Van-Leer (1977), "Towards the Ultimate Conservative Difference Scheme. IV. A New Approach to Numerical Convection"
613	!----------------------------------------------------------------------------------------------------------------------!
614	subroutine vl1d(nlayer,q,pente_max,masse,w,qm)
615
616	implicit none
617	!----------------------------------------------------------------------------------------------------------------------!
618	! VARIABLE DECLARATION
619	!----------------------------------------------------------------------------------------------------------------------!
620	! Input arguments:
621	!-----------------
622	integer, intent(in) :: &
623	nlayer ! number of layers
624
625	real, intent(in) :: &
626	pente_max, &! coefficient, pente_max = 2 advised
627	masse(nlayer), &! masse layer Dp/g (kg)
628	q(nlayer) ! quantity of tracer
629	!----------------------------------------------------------------------------------------------------------------------!
630	! In-Output arguments:
631	!---------------------
632	real, intent(inout) :: &
633	w(nlayer+1) ! masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
634	!----------------------------------------------------------------------------------------------------------------------!
635	! Output arguments:
636	!------------------
637	real, intent(out) :: &
638	qm(nlayer+1) ! quantity of tracer after Van-Leer scheme
639	!----------------------------------------------------------------------------------------------------------------------!
640	! Locals variables:
641	!------------------
642	integer :: &
643	l, &! loop on nlayer
644	m ! index
645
646	real :: &
647	dzqmax, &! maximum of dzq between two adjacent layers
648	sigw, &!
649	Mtot, &!
650	MQtot, &!
651	dzq(nlayer), &!
652	dzqw(nlayer),&!
653	adzqw(nlayer) !
654	!======================================================================================================================!
655	! BEGIN ===============================================================================================================!
656	!======================================================================================================================!
657	! 1. On oriente tout dans le sens de la pression: w > 0 when down
658	!----------------------------------------------------------------------------------------------------------------------!
659	do l = 2, nlayer
660	dzqw(l) = q(l-1) - q(l)
661	adzqw(l) = abs(dzqw(l))
662	end do
663
664	do l = 2, nlayer-1
665	if(dzqw(l)*dzqw(l+1)>0.) then
666	dzq(l) = 0.5 * (dzqw(l)+dzqw(l+1))
667	else
668	dzq(l) = 0.
669	end if
670
671	dzqmax = pente_max * min(adzqw(l), adzqw(l+1))
672
673	dzq(l) = sign(min(abs(dzq(l)),dzqmax), dzq(l))
674	end do
675
676	dzq(1)=0.
677	dzq(nlayer)=0.
678
679	do l = 1, nlayer-1
680	!----------------------------------------------------------------------------------------------------------------------!
681	! 2.1. Regular scheme (transfered mass < layer mass)
682	!----------------------------------------------------------------------------------------------------------------------!
683	if (w(l+1)>0. .and. w(l+1)<=masse(l+1)) then
684	sigw = w(l+1) / masse(l+1)
685	qm(l+1) = (q(l+1) + 0.5(1.-sigw)dzq(l+1))
686	else if (w(l+1)<=0. .and. -w(l+1)<=masse(l)) then
687	sigw = w(l+1) / masse(l)
688	qm(l+1) = (q(l) - 0.5(1.+sigw)dzq(l))
689	!----------------------------------------------------------------------------------------------------------------------!
690	! 2.2. Extended scheme (transfered mass > layer mass)
691	!----------------------------------------------------------------------------------------------------------------------!
692	else if (w(l+1)>0.) then
693	m = l+1
694	Mtot = masse(m)
695	MQtot = masse(m)*q(m)
696
697	do while ((m<nlayer).and.(w(l+1)>(Mtot+masse(m+1))))
698	m = m+1
699	Mtot = Mtot + masse(m)
700	MQtot = MQtot + masse(m)*q(m)
701	end do
702
703	if (m<nlayer) then
704	sigw = (w(l+1)-Mtot) / masse(m+1)
705	qm(l+1) = (1/w(l+1))( MQtot + (w(l+1)-Mtot) (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
706	else
707	w(l+1) = Mtot
708	qm(l+1) = Mqtot / Mtot
709	call abort_physic('co2condens4micro', 'top layer is disapearing !', 1)
710	end if
711	!----------------------------------------------------------------------------------------------------------------------!
712	else ! if(w(l+1).lt.0)
713	m = l-1
714	Mtot = masse(m+1)
715	MQtot = masse(m+1)*q(m+1)
716	if (m>0) then ! because some compilers will have problems evaluating masse(0)
717	do while ((m>0).and.(-w(l+1)>(Mtot+masse(m))))
718	m = m-1
719	Mtot = Mtot + masse(m+1)
720	MQtot = MQtot + masse(m+1)*q(m+1)
721	if (m==0) then
722	exit
723	end if
724	end do
725	end if
726
727	if (m>0) then
728	sigw = (w(l+1)+Mtot) / masse(m)
729	qm(l+1) = (-1/w(l+1)) * ( MQtot + (-w(l+1)-Mtot) * (q(m)-0.5(1.+sigw)dzq(m)) )
730	else
731	qm(l+1) = (-1/w(l+1)) * (MQtot + (-w(l+1)-Mtot)*qm(1))
732	end if
733	end if
734	end do ! l = 1, nlayer-1
735	!======================================================================================================================!
736	! END =================================================================================================================!
737	!======================================================================================================================!
738	end subroutine vl1d
739
740	end module co2condens_mod4micro

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

Download in other formats:

Original Format