Context Navigation

newcondens.F @ 1818

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

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

Download in other formats:

Original Format