Context Navigation

newcondens.F @ 1009

Last change on this file since 1009 was 890, checked in by emillour, 12 years ago
Mars GCM: minor bug correction in newcondens.F and some code tidying (in co2snow.F also). EM
File size: 31.9 KB

Line
1	SUBROUTINE newcondens(ngrid,nlayer,nq,ptimestep,
2	$ pcapcal,pplay,pplev,ptsrf,pt,
3	$ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
4	$ piceco2,psolaralb,pemisurf,
5	$ pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
6	$ fluxsurf_sw,zls)
7
8	IMPLICIT NONE
9	c=======================================================================
10	c subject:
11	c --------
12	c Condensation/sublimation of CO2 ice on the ground and in the
13	c atmosphere
14	c (Scheme described in Forget et al., Icarus, 1998)
15	c
16	c author: Francois Forget 1994-1996
17	c ------
18	c
19	c input:
20	c -----
21	c ngrid nombre de points de verticales
22	c (toutes les boucles de la physique sont au
23	c moins vectorisees sur ngrid)
24	c nlayer nombre de couches
25	c pplay(ngrid,nlayer) Pressure levels
26	c pplev(ngrid,nlayer+1) Pressure levels
27	c pt(ngrid,nlayer) temperature (en K)
28	c ptsrf(ngrid) temperature de surface
29	c
30	c \
31	c pdt(ngrid,nlayermx) \ derivee temporelle physique avant condensation
32	c / ou sublimation pour pt,ptsrf
33	c pdtsrf(ngrid) /
34	c
35	c output:
36	c -------
37	c
38	c pdpsrf(ngrid) \ derivee temporelle physique (contribution de
39	c pdtc(ngrid,nlayermx) / la condensation ou sublimation) pour Ps,pt,ptsrf
40	c pdtsrfc(ngrid) /
41	c
42	c Entree/sortie
43	c -------------
44	c
45	c piceco2(ngrid) : quantite de glace co2 au sol (kg/m2)
46	c psolaralb(ngrid,2) : albedo au sol
47	c pemisurf(ngrid) : emissivite du sol
48
49	c
50	c=======================================================================
51	c
52	c 0. Declarations :
53	c ------------------
54	c
55	#include "dimensions.h"
56	#include "dimphys.h"
57	#include "comcstfi.h"
58	#include "surfdat.h"
59	#include "comgeomfi.h"
60	#include "comvert.h"
61	#include "paramet.h"
62	#include "callkeys.h"
63	#include "tracer.h"
64
65	c-----------------------------------------------------------------------
66	c Arguments :
67	c ---------
68	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
69	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
70	INTEGER,INTENT(IN) :: nq ! number of tracers
71
72	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
73	REAL,INTENT(IN) :: pcapcal(ngrid)
74	REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
75	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
76	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
77	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
78	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
79	REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
80	! previous physical processes (K/s)
81	REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
82	! from previous physical processes
83	REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
84	! from previous physical processes
85	REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from
86	! previous physical processes (K/s)
87	REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
88	REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
89	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
90	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
91	! previous physical processes
92	REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2)
93	REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface
94	REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface
95
96	! tendencies due to CO2 condensation/sublimation:
97	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
98	REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s)
99	REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
100	REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
101	REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
102	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
103
104	! added to calculate flux dependent albedo:
105	REAL,intent(in) :: fluxsurf_sw(ngrid,2)
106	real,intent(in) :: zls ! solar longitude (rad)
107
108	c
109	c Local variables :
110	c -----------------
111
112	c variables used for albedo parametrization
113	c --------------------------------------------
114	INTEGER i,j
115	REAL Fluxmean(jjp1)
116	INTEGER l,ig,iq,icap,nmix
117	LOGICAL transparency, fluxdependent
118	c flag transparency if you want to make the co2ice semi-transparent
119	PARAMETER(transparency=.true.)
120	c flag fluxdependent if you want the co2ice albedo to be dependent on
121	c the incident solar flux
122	PARAMETER(fluxdependent=.false.)
123	REAL slopy,alpha,constA,constB,constT,albediceF_new(ngrid)
124	REAL zt(ngrid,nlayer)
125	REAL zcpi
126	REAL ztcond (ngrid,nlayer+1)
127	REAL ztcondsol(ngrid)
128	REAL zdiceco2(ngrid)
129	REAL zcondicea(ngrid,nlayer)
130	REAL zcondices(ngrid)
131	REAL zfallice(ngrid,nlayer+1) , zfallheat
132	REAL zmflux(nlayer+1)
133	REAL zu(nlayer),zv(nlayer)
134	REAL zq(nlayer,nq),zq1(nlayer)
135	REAL ztsrf(ngrid)
136	REAL ztc(nlayer), ztm(nlayer+1)
137	REAL zum(nlayer+1) , zvm(nlayer+1)
138	REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
139	REAL masse(nlayer),w(nlayer+1)
140	REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
141	LOGICAL condsub(ngrid)
142
143	c variable speciale diagnostique
144	real tconda1(ngrid,nlayer)
145	real tconda2(ngrid,nlayer)
146	c REAL zdiceco2a(ngrid) ! for diagnostic only
147	real zdtsig (ngrid,nlayer)
148	real zdt (ngrid,nlayer)
149	real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
150	! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
151	! then condensation temperature is computed using partial pressure of CO2
152	logical,parameter :: improved_ztcond=.true.
153	! Bound co2 (tracer) values...
154	logical,parameter :: bound_qco2=.false.
155	real,parameter :: qco2max=1.1
156	real,parameter :: qco2mini=0.1
157	real :: zqco2
158
159	c local saved variables
160	integer,save :: ico2 ! index of CO2 tracer
161	real,save :: qco2min,qco2,mmean
162	real,allocatable,save :: emisref(:)
163	real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
164	real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
165	real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
166	REAL,SAVE :: acond,bcond,ccond
167	! REAL,SAVE :: albediceF(ngridmx)
168	real,save :: m_co2, m_noco2, A , B
169
170	LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
171
172	integer flag
173
174	c----------------------------------------------------------------------
175
176	c Initialisation
177	c --------------
178	c
179	IF (firstcall) THEN
180
181	allocate(emisref(ngrid))
182
183	bcond=1./tcond1mb
184	ccond=cpp/(g*latcond)
185	acond=r/latcond
186
187	firstcall=.false.
188	write(,) 'Newcondens: improved_ztcond=',improved_ztcond
189	write(,) 'Newcondens: bound_qco2=',bound_qco2
190	PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
191	PRINT*,'acond,bcond,ccond',acond,bcond,ccond
192
193	ico2=0
194
195	if (tracer) then
196	c Prepare Special treatment if one of the tracer is CO2 gas
197	do iq=1,nqmx
198	if (noms(iq).eq."co2") then
199	ico2=iq
200	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
201	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
202	c Compute A and B coefficient use to compute
203	c mean molecular mass Mair defined by
204	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
205	c 1/Mair = A*q(ico2) + B
206	A =(1/m_co2 - 1/m_noco2)
207	B=1/m_noco2
208	endif
209	enddo
210	c minimum CO2 mix. ratio below which mixing occur with layer above:
211	qco2min =0.75
212	end if
213	ENDIF ! of IF (firstcall)
214	zcpi=1./cpp
215	c
216	c======================================================================
217	c Calcul of CO2 condensation sublimation
218	c ============================================================
219	c
220	c Used variable :
221	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
222	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
223	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
224	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
225	c
226	c pdtc(ngrid,nlayermx) : dT/dt due to cond/sub
227	c
228	c
229	c Tendencies set to 0 (except pdtc)
230	c -------------------------------------
231	DO l=1,nlayer
232	DO ig=1,ngrid
233	zcondicea(ig,l) = 0.
234	zfallice(ig,l) = 0.
235	pduc(ig,l) = 0
236	pdvc(ig,l) = 0
237	END DO
238	END DO
239
240	DO iq=1,nqmx
241	DO l=1,nlayer
242	DO ig=1,ngrid
243	pdqc(ig,l,iq) = 0
244	END DO
245	END DO
246	END DO
247
248	DO ig=1,ngrid
249	zfallice(ig,nlayer+1) = 0.
250	zcondices(ig) = 0.
251	pdtsrfc(ig) = 0.
252	pdpsrf(ig) = 0.
253	condsub(ig) = .false.
254	zdiceco2(ig) = 0.
255	ENDDO
256	zfallheat=0
257
258	c *************************
259	c ATMOSPHERIC CONDENSATION
260	c *************************
261
262	c Compute CO2 Volume mixing ratio
263	c -------------------------------
264	if (improved_ztcond.and.(ico2.ne.0)) then
265	DO l=1,nlayer
266	DO ig=1,ngrid
267	qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
268	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
269	mmean=1/(A*qco2 +B)
270	vmr_co2(ig,l) = qco2*mmean/m_co2
271	ENDDO
272	ENDDO
273	else
274	DO l=1,nlayer
275	DO ig=1,ngrid
276	vmr_co2(ig,l)=0.95
277	ENDDO
278	ENDDO
279	end if
280
281	c forecast of atmospheric temperature zt and frost temperature ztcond
282	c --------------------------------------------------------------------
283
284	DO l=1,nlayer
285	DO ig=1,ngrid
286	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
287	! ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
288	ztcond(ig,l)=
289	& 1./(bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
290	if (pplay(ig,l).lt.1e-4) ztcond(ig,l)=0.0 !mars Monica
291	ENDDO
292	ENDDO
293
294	ztcond(:,nlayer+1)=ztcond(:,nlayer)
295
296	c Condensation/sublimation in the atmosphere
297	c ------------------------------------------
298	c (calcul of zcondicea , zfallice and pdtc)
299	c
300	DO l=nlayer , 1, -1
301	DO ig=1,ngrid
302	pdtc(ig,l)=0.
303	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
304	condsub(ig)=.true.
305	IF (zfallice(ig,l+1).gt.0) then
306	zfallheat=zfallice(ig,l+1)*
307	& (pphi(ig,l+1)-pphi(ig,l) +
308	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
309	ELSE
310	zfallheat=0.
311	ENDIF
312	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
313	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
314	& ccondpdtc(ig,l)- zfallheat
315	c Case when the ice from above sublimes entirely
316	c """""""""""""""""""""""""""""""""""""""""""""""
317	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
318	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
319	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
320	zcondicea(ig,l)= -zfallice(ig,l+1)
321	END IF
322
323	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
324	END IF
325	ENDDO
326	ENDDO
327
328	c *************************
329	c SURFACE CONDENSATION
330	c *************************
331
332	c forecast of ground temperature ztsrf and frost temperature ztcondsol
333	c --------------------------------------------------------------------
334	DO ig=1,ngrid
335	ztcondsol(ig)=
336	& 1./(bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
337	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
338	ENDDO
339
340	c
341	c Condensation/sublimation on the ground
342	c --------------------------------------
343	c (calcul of zcondices , pdtsrfc)
344	c
345	DO ig=1,ngrid
346	IF(ig.GT.ngrid/2+1) THEN
347	icap=2
348	ELSE
349	icap=1
350	ENDIF
351
352	c
353	c Loop on where we have condensation/ sublimation
354	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
355	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
356	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
357	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
358	condsub(ig) = .true.
359
360	IF (zfallice(ig,1).gt.0) then
361	zfallheat=zfallice(ig,1)*
362	& (pphi(ig,1)- phisfi(ig) +
363	& cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
364	ELSE
365	zfallheat=0.
366	ENDIF
367
368	c condensation or partial sublimation of CO2 ice
369	c """""""""""""""""""""""""""""""""""""""""""""""
370	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
371	& /(latcond*ptimestep) - zfallheat
372	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
373
374	c If the entire CO_2 ice layer sublimes
375	c """"""""""""""""""""""""""""""""""""""""""""""""""""
376	c (including what has just condensed in the atmosphere)
377
378	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
379	& -zcondices(ig))THEN
380	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
381	pdtsrfc(ig)=(latcond/pcapcal(ig))*
382	& (zcondices(ig)+zfallheat)
383	END IF
384
385	c Changing CO2 ice amount and pressure :
386	c """"""""""""""""""""""""""""""""""""
387
388	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
389	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
390	pdpsrf(ig) = -zdiceco2(ig)*g
391
392	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
393	PRINT*,'STOP in condens'
394	PRINT*,'condensing more than total mass'
395	PRINT*,'Grid point ',ig
396	PRINT*,'Ps = ',pplev(ig,1)
397	PRINT*,'d Ps = ',pdpsrf(ig)
398	STOP
399	ENDIF
400	END IF ! if there is condensation/sublimmation
401	ENDDO ! of DO ig=1,ngrid
402
403	c ********************************************************************
404	c Surface albedo and emissivity of the surface below the snow (emisref)
405	c ********************************************************************
406	c Prepare the case where albedo varies with insolation:
407	c ----------------------------------------------------
408	! if (fluxdependent) then
409	!
410	c Calcul du flux moyen (zonal mean)
411	! do j=1,jjp1
412	! Fluxmean(j)=0
413	! do i=1,iim
414	! ig=1+(j-2)*iim +i
415	! if(j.eq.1) ig=1
416	! if(j.eq.jjp1) ig=ngrid
417	! Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
418	! $ +fluxsurf_sw(ig,2)
419	! enddo
420	! Fluxmean(j)=Fluxmean(j)/float(iim)
421	! enddo
422	!
423	c const A and B used to calculate the albedo which depends on solar flux
424	c albedice=constA+constB*Flux
425	c constT = time step to calculate the solar flux when flux decreases
426	! constA=0.26
427	c constA=0.33
428	c constA=0.186
429	! constB=0.00187
430	! constT=10
431	! endif ! of if (fluxdependent)
432
433	! Check that amont of CO2 ice is not problematic
434	DO ig=1,ngrid
435	if(.not.piceco2(ig).ge.0.) THEN
436	if(piceco2(ig).le.-5.e-8) print*,
437	$ 'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
438	piceco2(ig)=0.
439	endif
440	ENDDO
441
442	! Set albedo and emissivity of the surface
443	! ----------------------------------------
444	CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
445
446	c Calcul de l'albedo
447	c ------------------
448	! do ig =1,ngrid
449	! IF(ig.GT.ngrid/2+1) THEN
450	! icap=2
451	! ELSE
452	! icap=1
453	! ENDIF
454	! IF(firstcall2) THEN
455	! albediceF(ig)=albedice(icap)
456	! ENDIF
457	c if there is still co2ice ccccccccccccccccccccccc
458	! if (piceco2(ig).gt.0) then
459	! emisref(ig) = emisice(icap)
460
461	c if flux dependent albedo is used
462	c --------------------------------
463	! if (fluxdependent) then
464	! j=INT((ig-2)/iim)+2
465	! if(ig.eq.1) j=1
466	! if(ig.eq.ngrid) j=jjp1
467	c albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
468	c $ constA+constB*250)
469	! albediceF_new(ig)=constA+constB*Fluxmean(j)
470	! if (albediceF(ig).gt.albediceF_new(ig)) then
471	! albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
472	! $ constT)*(albediceF_new(ig)-albediceF(ig))
473	! else
474	! albediceF(ig)=albediceF_new(ig)
475	! endif
476	c if part of the ice is transparent
477	c slopy = pente de la droite: alpha = y*co2ice/1620
478	c pour une valeur superieur a une epaisseur de glace donnee
479	c ici, epaisseur limite = 10cm
480	! if (transparency) then
481	! slopy=1/(1620*0.10)
482	! alpha=MIN(slopy*piceco2(ig),1.)
483	! psolaralb(ig,1) = alpha*albediceF(ig)
484	! $ +(1-alpha)*albedodat(ig)
485	! psolaralb(ig,2) = psolaralb(ig,1)
486	! else
487	! psolaralb(ig,1) = albediceF(ig)
488	! psolaralb(ig,2) = psolaralb(ig,1)
489	! endif
490	! else
491	c transparency set to true and fluxdependent set to false
492	! if (transparency) then
493	! slopy=1/(1620*0.10)
494	! alpha=MIN(slopy*piceco2(ig),1.)
495	! psolaralb(ig,1) = alpha*albedice(icap)
496	! $ +(1-alpha)*albedodat(ig)
497	! psolaralb(ig,2) = psolaralb(ig,1)
498	! else
499	c simplest case: transparency and flux dependent set to false
500	! psolaralb(ig,1) = albedice(icap)
501	! psolaralb(ig,2) = albedice(icap)
502	! endif
503	! endif
504	c no more co2ice, albedo = ground albedo
505	! else
506	! psolaralb(ig,1) = albedodat(ig)
507	! psolaralb(ig,2) = albedodat(ig)
508	! emisref(ig) = emissiv
509	! pemisurf(ig) = emissiv
510	! endif
511	! end do ! end of the ig loop
512
513	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
514	DO ig=1,ngrid
515	if (piceco2(ig).eq.0) then
516	pemisurf(ig)=emissiv
517	endif
518	ENDDO
519
520	! firstcall2=.false.
521	c ***************************************************************
522	c Correction to account for redistribution between sigma or hybrid
523	c layers when changing surface pressure (and warming/cooling
524	c of the CO2 currently changing phase).
525	c *************************************************************
526
527	DO ig=1,ngrid
528	if (condsub(ig)) then
529	do l=1,nlayer
530	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
531	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
532	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
533	do iq=1,nqmx
534	zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
535	enddo
536	end do
537
538	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
539	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
540
541	zmflux(1) = -zcondices(ig)
542	DO l=1,nlayer
543	zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
544	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
545	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
546	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
547	END DO
548
549	c Mass of each layer
550	c ------------------
551	DO l=1,nlayer
552	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
553	END DO
554
555
556	c Corresponding fluxes for T,U,V,Q
557	c """"""""""""""""""""""""""""""""
558
559	c averaging operator for TRANSPORT
560	c """"""""""""""""""""""""""""""""
561	c Value transfert at the surface interface when condensation
562	c sublimation:
563	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
564	zum(1) = 0
565	zvm(1) = 0
566	do iq=1,nqmx
567	zqm(1,iq)=0. ! most tracer do not condense !
568	enddo
569	c Special case if one of the tracer is CO2 gas
570	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
571
572	c Van Leer scheme:
573	DO l=1,nlayer+1
574	w(l)=-zmflux(l)*ptimestep
575	END DO
576	call vl1d(ztc,2.,masse,w,ztm)
577	call vl1d(zu ,2.,masse,w,zum)
578	call vl1d(zv ,2.,masse,w,zvm)
579	do iq=1,nqmx
580	do l=1,nlayer
581	zq1(l)=zq(l,iq)
582	enddo
583	zqm1(1)=zqm(1,iq)
584	call vl1d(zq1,2.,masse,w,zqm1)
585	do l=2,nlayer
586	zq( l,iq)=zq1(l)
587	zqm(l,iq)=zqm1(l)
588	enddo
589	enddo
590
591	c Surface condensation affects low winds
592	if (zmflux(1).lt.0) then
593	zum(1)= zu(1) * (w(1)/masse(1))
594	zvm(1)= zv(1) * (w(1)/masse(1))
595	if (w(1).gt.masse(1)) then ! ensure numerical stability
596	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
597	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
598	end if
599	end if
600
601	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
602	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
603	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
604	do iq=1,nqmx
605	zqm(nlayer+1,iq)= zq(nlayer,iq)
606	enddo
607
608	#ifdef MESOSCALE
609	!!!! AS: This part must be commented in the mesoscale model
610	!!!! AS: ... to avoid instabilities.
611	!!!! AS: you have to compile with -DMESOSCALE to do so
612	#else
613	c Tendencies on T, U, V, Q
614	c """"""""""""""""""""""""
615	DO l=1,nlayer
616
617	c Tendencies on T
618	zdtsig(ig,l) = (1/masse(l)) *
619	& ( zmflux(l)*(ztm(l) - ztc(l))
620	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
621	& + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
622	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
623
624	c Tendencies on U
625	pduc(ig,l) = (1/masse(l)) *
626	& ( zmflux(l)*(zum(l) - zu(l))
627	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
628
629
630	c Tendencies on V
631	pdvc(ig,l) = (1/masse(l)) *
632	& ( zmflux(l)*(zvm(l) - zv(l))
633	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
634
635	END DO
636	#endif
637
638	c Tendencies on Q
639	do iq=1,nqmx
640	! if (noms(iq).eq.'co2') then
641	if (iq.eq.ico2) then
642	c SPECIAL Case when the tracer IS CO2 :
643	DO l=1,nlayer
644	pdqc(ig,l,iq)= (1/masse(l)) *
645	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
646	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
647	& + zcondicea(ig,l)*(zq(l,iq)-1.) )
648	END DO
649	else
650	DO l=1,nlayer
651	pdqc(ig,l,iq)= (1/masse(l)) *
652	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
653	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
654	& + zcondicea(ig,l)*zq(l,iq) )
655	END DO
656	end if
657	enddo
658
659	c --------------------------------------------------------
660	c Roughly Simulate Molecular mixing when CO2 is too depleted by
661	c Surface condensation (mixing starts if qco2 < qco2min )
662	c FF 06/2004
663	c WARNING : this is now done in convadj, better (FF 02/2005)
664	c --------------------------------------------------------
665	flag=0 ! now done in convadj : must be =0
666	if (flag.eq.1) then
667	if(ico2.gt.0) then ! relevant only if one tracer is CO2
668	if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
669	& .lt.qco2min) then
670	do iq=1,nqmx
671	zq(1,iq)=pq(ig,1,iq)
672	& +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
673	Smq(1,iq) = masse(1)*zq(1,iq)
674	end do
675	Sm(1) = masse(1)
676	do l =2,nlayermx
677	do iq=1,nqmx
678	zq(l,iq)=pq(ig,l,iq)
679	& +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
680	smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
681	end do
682	sm(l) = sm(l-1) + masse(l)
683	if(zq(l,ico2).gt.qco2min) then
684	c mixmas: mass of atmosphere that must be mixed to reach qco2min
685	mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
686	& /(zq(l,ico2)-qco2min)
687	if((mixmas.le.sm(l)))then
688	c OK if mixed mass less than mass of layers affected
689	nmix=l ! number of layer affected by mixing
690	goto 99
691	end if
692	end if
693	end do
694	99 continue
695	do iq=1,nqmx
696	qmix=zq(nmix,iq)
697	& +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
698	do l=1,nmix-1
699	pdqc(ig,l,iq)=
700	& (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
701	end do
702	c layer only partly mixed :
703	pdqc(ig,nmix,iq)=(
704	& qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
705	& -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
706
707	end do
708	end if
709	end if
710
711	endif ! (flag.eq.1)
712	end if ! if (condsub)
713	END DO ! loop on ig
714
715	c ***************************************************************
716	c CO2 snow / clouds scheme
717	c ***************************************************************
718
719	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
720	& zcondicea,zcondices,zfallice,pemisurf)
721
722	c ***************************************************************
723	c Ecriture des diagnostiques
724	c ***************************************************************
725
726	c DO l=1,nlayer
727	c DO ig=1,ngrid
728	c Taux de cond en kg.m-2.pa-1.s-1
729	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
730	c Taux de cond en kg.m-3.s-1
731	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
732	c END DO
733	c END DO
734	c call WRITEDIAGFI(ngridmx,'tconda1',
735	c &'Taux de condensation CO2 atmospherique /Pa',
736	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
737	c call WRITEDIAGFI(ngridmx,'tconda2',
738	c &'Taux de condensation CO2 atmospherique /m',
739	c & 'kg.m-3.s-1',3,tconda2)
740
741	! output falling co2 ice in 1st layer:
742	! call WRITEDIAGFI(ngridmx,'fallice',
743	! &'Precipitation of co2 ice',
744	! & 'kg.m-2.s-1',2,zfallice(1,1))
745
746	!! Specific stuff to bound co2 tracer values ....
747	if (bound_qco2.and.(ico2.ne.0)) then
748	do ig=1,ngrid
749	do l=1,nlayer
750	zqco2=pq(ig,l,ico2)
751	& +(pdq(ig,l,ico2)+pdqc(ig,l,ico2))*ptimestep
752	if (zqco2.gt.qco2max) then
753	! correct pdqc:
754	pdqc(ig,l,ico2)=((qco2max-pq(ig,l,ico2))/ptimestep)
755	& -pdq(ig,l,ico2)
756	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
757	& " so that co2 conc. does not exceed",qco2max
758	write(,) " ig:",ig," l:",l
759	endif ! of if (zqco2.gt.qco2max)
760	if (zqco2.lt.qco2mini) then
761	! correct pdqc:
762	pdqc(ig,l,ico2)=((qco2mini-pq(ig,l,ico2))/ptimestep)
763	& -pdq(ig,l,ico2)
764	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
765	& " so that co2 conc. is not less than",qco2mini
766	write(,) " ig:",ig," l:",l
767	endif ! of if (zqco2.lt.qco2mini)
768	end do
769	enddo
770	endif ! of if (bound_qco2.and.(ico2.ne.0)) then
771
772	return
773	end
774
775
776
777	c *****************************************************************
778	SUBROUTINE vl1d(q,pente_max,masse,w,qm)
779	c
780	c
781	c Operateur de moyenne inter-couche pour calcul de transport type
782	c Van-Leer " pseudo amont " dans la verticale
783	c q,w sont des arguments d'entree pour le s-pg ....
784	c masse : masse de la couche Dp/g
785	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
786	c pente_max = 2 conseillee
787	c
788	c
789	c --------------------------------------------------------------------
790	IMPLICIT NONE
791
792	#include "dimensions.h"
793
794	c
795	c
796	c
797	c Arguments:
798	c ----------
799	real masse(llm),pente_max
800	REAL q(llm),qm(llm+1)
801	REAL w(llm+1)
802	c
803	c Local
804	c ---------
805	c
806	INTEGER l
807	c
808	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
809	real sigw, Mtot, MQtot
810	integer m
811	c integer ismax,ismin
812
813
814	c On oriente tout dans le sens de la pression
815	c W > 0 WHEN DOWN !!!!!!!!!!!!!
816
817	do l=2,llm
818	dzqw(l)=q(l-1)-q(l)
819	adzqw(l)=abs(dzqw(l))
820	enddo
821
822	do l=2,llm-1
823	if(dzqw(l)*dzqw(l+1).gt.0.) then
824	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
825	else
826	dzq(l)=0.
827	endif
828	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
829	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
830	enddo
831
832	dzq(1)=0.
833	dzq(llm)=0.
834
835	do l = 1,llm-1
836
837	c Regular scheme (transfered mass < layer mass)
838	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
839	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
840	sigw=w(l+1)/masse(l+1)
841	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
842	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
843	sigw=w(l+1)/masse(l)
844	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
845
846	c Extended scheme (transfered mass > layer mass)
847	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
848	else if(w(l+1).gt.0.) then
849	m=l+1
850	Mtot = masse(m)
851	MQtot = masse(m)*q(m)
852	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1))))
853	m=m+1
854	Mtot = Mtot + masse(m)
855	MQtot = MQtot + masse(m)*q(m)
856	end do
857	if (m.lt.llm) then
858	sigw=(w(l+1)-Mtot)/masse(m+1)
859	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
860	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
861	else
862	w(l+1) = Mtot
863	qm(l+1) = Mqtot / Mtot
864	write(,) 'top layer is disapearing !'
865	stop
866	end if
867	else ! if(w(l+1).lt.0)
868	m = l-1
869	Mtot = masse(m+1)
870	MQtot = masse(m+1)*q(m+1)
871	if (m.gt.0) then ! because some compilers will have problems
872	! evaluating masse(0)
873	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
874	m=m-1
875	Mtot = Mtot + masse(m+1)
876	MQtot = MQtot + masse(m+1)*q(m+1)
877	if (m.eq.0) exit
878	end do
879	endif
880	if (m.gt.0) then
881	sigw=(w(l+1)+Mtot)/masse(m)
882	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
883	& (q(m)-0.5(1.+sigw)dzq(m)) )
884	else
885	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
886	end if
887	end if
888	enddo
889
890	c boundary conditions (not used in newcondens !!)
891	c qm(llm+1)=0.
892	c if(w(1).gt.0.) then
893	c qm(1)=q(1)
894	c else
895	c qm(1)=0.
896	c end if
897
898	return
899	end

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

Download in other formats:

Original Format