Context Navigation

newcondens.F @ 86

Last change on this file since 86 was 86, checked in by aslmd, 14 years ago

*
mars + LMD_MM_MARS * Precompilation flag MESOSCALE for better transparency

* in shared phymars between GCM and mesoscale model

M 85 mars/libf/phymars/meso_physiq.F
M 85 mars/libf/phymars/meso_inifis.F
Added a pre-compilation flag MESOSCALE so that the LMDZ.MARS GCM
will compile without stating errors because of mesoscale routines.

M 85 mars/libf/phymars/newcondens.F
M 85 mars/libf/phymars/testphys1d.F
M 85 mars/libf/phymars/dustlift.F
D 85 mars/libf/phymars/meso_testphys1d.F
D 85 mars/libf/phymars/meso_dustlift.F
D 85 mars/libf/phymars/meso_newcondens.F
Now, this MESOSCALE precompilation flag can be used to lower
the number of meso_* routines when adaptations for mesoscale
applications are not very extended.
--> Three meso_* routines were deleted and changes are
now impacted under the MESOSCALE flag in the original GCM routines
--> Completely transparent for GCM compilation since it is devoid
of the -DMESOSCALE option
--> Very good for syncing because changes in dustlift, newcondens
will be directly available in the mesoscale model

M 84 mesoscale/LMD_MM_MARS/makemeso
Changed meso_testphys1d in testphys1d

M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_pgf
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpifort
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_ifort
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_g95
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpi
Added the option -DMESOSCALE in these scripts

*
LMD_MM_MARS * Various minor changes related to water cycle and plotting routines

* Also included the GW test case

A 0 mesoscale/LMDZ.MARS.new/myGCM/DEFS_JB/callphys.def.orig
M 84 mesoscale/NOTES.txt
D 84 mesoscale/LMD_MM_MARS/SRC/ARWpost/idl
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
M 84 mesoscale/LMD_MM_MARS/SIMU/gnome_launch.meso
M 85 mesoscale/PLOT/MINIMAL/map_latlon.pro
D 85 mesoscale/PLOT/SPEC/LES/getget.pro
M 85 mesoscale/PLOT/SPEC/MAP/map_uvt.pro
A + - mesoscale/PLOT/SPEC/getget.pro
A 0 mesoscale/PLOT/RESERVE/obsolete
A 0 mesoscale/TESTS/TESTGW.tar.gz
M 84 000-USERS

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

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

Download in other formats:

Original Format