Context Navigation

newcondens.F @ 1292

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

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

Download in other formats:

Original Format