Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

condens_n2.F90 @ 3541

Last change on this file since 3541 was 3175, checked in by emillour, 18 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 30.9 KB

Line
1	subroutine condens_n2(ngrid,nlayer,nq,ptimestep, &
2	pcapcal,pplay,pplev,ptsrf,pt, &
3	pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
4	picen2,psolaralb,pemisurf, &
5	pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
6	pdqc,pdicen2)
7
8	use radinc_h, only : naerkind
9	use comgeomfi_h
10	implicit none
11
12	!==================================================================
13	! Purpose
14	! -------
15	! Condense and/or sublime N2 ice on the ground and in the
16	! atmosphere, and sediment the ice.
17	!
18	! Inputs
19	! ------
20	! ngrid Number of vertical columns
21	! nlayer Number of layers
22	! pplay(ngrid,nlayer) Pressure layers
23	! pplev(ngrid,nlayer+1) Pressure levels
24	! pt(ngrid,nlayer) Temperature (in K)
25	! ptsrf(ngrid) Surface temperature
26	!
27	! pdt(ngrid,nlayermx) Time derivative before condensation/sublimation of pt
28	! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf
29	! picen2(ngrid) n2 ice at the surface (kg/m2)
30	!
31	! Outputs
32	! -------
33	! pdpsrf(ngrid) \ Contribution of condensation/sublimation
34	! pdtc(ngrid,nlayermx) / to the time derivatives of Ps, pt, and ptsrf
35	! pdtsrfc(ngrid) /
36	! pdicen2(ngrid) Tendancy n2 ice at the surface (kg/m2)
37	!
38	! Both
39	! ----
40	! psolaralb(ngrid) Albedo at the surface
41	! pemisurf(ngrid) Emissivity of the surface
42	!
43	! Authors
44	! -------
45	! Francois Forget (1996,2013)
46	! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
47	!
48	!
49	!==================================================================
50
51	#include "dimensions.h"
52	#include "dimphys.h"
53	#include "comcstfi.h"
54	#include "surfdat.h"
55	#include "comvert.h"
56	#include "callkeys.h"
57	#include "tracer.h"
58
59	!-----------------------------------------------------------------------
60	! Arguments
61
62	INTEGER ngrid, nlayer, nq
63
64	REAL ptimestep
65	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
66	REAL pcapcal(ngrid)
67	REAL pt(ngrid,nlayer)
68	REAL ptsrf(ngrid),flu1(ngrid),flu2(ngrid),flu3(ngrid)
69	REAL pphi(ngrid,nlayer)
70	REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
71	REAL pdtsrfc(ngrid),pdpsrf(ngrid)
72	REAL picen2(ngrid),psolaralb(ngrid),pemisurf(ngrid)
73
74
75	REAL pu(ngrid,nlayer) , pv(ngrid,nlayer)
76	REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
77	REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
78	REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
79	REAL pdqc(ngrid,nlayer,nq)
80
81	!-----------------------------------------------------------------------
82	! Local variables
83
84	INTEGER l,ig,ilay,it,iq,ich4_gas
85
86	REAL*8 zt(ngridmx,nlayermx)
87	REAL tcond_n2
88	REAL pcond_n2
89	REAL glob_average2d ! function
90	REAL zqn2(ngridmx,nlayermx) ! N2 MMR used to compute Tcond/zqn2
91	REAL ztcond (ngridmx,nlayermx)
92	REAL ztcondsol(ngridmx),zfallheat
93	REAL pdicen2(ngridmx)
94	REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
95	REAL zfallice(ngridmx,nlayermx+1)
96	REAL zmflux(nlayermx+1)
97	REAL zu(nlayermx),zv(nlayermx)
98	REAL zq(nlayermx,nqmx),zq1(nlayermx)
99	REAL ztsrf(ngridmx)
100	REAL ztc(nlayermx), ztm(nlayermx+1)
101	REAL zum(nlayermx+1) , zvm(nlayermx+1)
102	REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
103	LOGICAL condsub(ngridmx)
104	REAL subptimestep
105	Integer Ntime
106	real masse (nlayermx),w(nlayermx+1)
107	real wq(ngridmx,nlayermx+1)
108	real vstokes,reff
109	real dWtotsn2
110	real condnconsn2(ngridmx)
111	real nconsMAXn2
112	! Special diagnostic variables
113	real tconda1(ngridmx,nlayermx)
114	real tconda2(ngridmx,nlayermx)
115	real zdtsig (ngridmx,nlayermx)
116	real zdtlatent (ngridmx,nlayermx)
117	real zdt (ngridmx,nlayermx)
118	REAL albediceF(ngridmx)
119	SAVE albediceF
120	INTEGER nsubtimestep,itsub !number of subtimestep when calling vl1d
121	REAL subtimestep !ptimestep/nsubtimestep
122	REAL dtmax
123
124	REAL zplevhist(ngridmx)
125	REAL zplevnew(ngridmx)
126	REAL zplev(ngridmx)
127	REAL zpicen2(ngridmx)
128	REAL ztsrfhist(ngridmx)
129	REAL zdtsrf(ngridmx)
130	REAL globzplevnew
131
132	REAL vmrn2(ngridmx)
133	SAVE vmrn2
134	REAL stephan
135	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
136	SAVE stephan
137	!-----------------------------------------------------------------------
138	! Saved local variables
139
140	! REAL latcond
141	REAL ccond
142	REAL cpice ! for atm condensation
143	SAVE cpice
144	! SAVE latcond,ccond
145	SAVE ccond
146
147	LOGICAL firstcall
148	SAVE firstcall
149	REAL SSUM
150	EXTERNAL SSUM
151
152	! DATA latcond /2.5e5/
153	! DATA latcond /1.98e5/
154	DATA cpice /1300./
155	DATA firstcall/.true./
156
157	INTEGER,SAVE :: i_n2ice=0 ! n2 ice
158	CHARACTER(LEN=20) :: tracername ! to temporarily store text
159	logical olkin ! option to prevent N2 ice effect in the south
160	DATA olkin/.false./
161	save olkin
162
163	CHARACTER(len=10) :: tname
164
165	!-----------------------------------------------------------------------
166
167	! Initialisation
168	IF (firstcall) THEN
169	ccond=cpp/(g*lw_n2)
170	print*,'In condens_n2cloud: ccond=',ccond,' latcond=',lw_n2
171
172	! calculate global mean surface pressure for the fast mode
173	IF (fast) THEN
174	DO ig=1,ngridmx
175	kp(ig) = exp(-phisfi(ig)/(r*38.))
176	ENDDO
177	p00=glob_average2d(kp) ! mean pres at ref level
178	ENDIF
179
180	vmrn2(:) = 1.
181	IF (ch4lag) then
182	DO ig=1,ngridmx
183	if (lati(ig)*180./pi.ge.latlag) then
184	vmrn2(ig) = vmrlag
185	endif
186	ENDDO
187	ENDIF
188	IF (no_n2frost) then
189	DO ig=1,ngridmx
190	if (picen2(ig).eq.0.) then
191	vmrn2(ig) = 1.e-15
192	endif
193	ENDDO
194	ENDIF
195	firstcall=.false.
196	ENDIF
197
198	!-----------------------------------------------------------------------
199	! Calculation of n2 condensation / sublimation
200
201	! Variables used:
202	! picen2(ngrid) : amount of n2 ice on the ground (kg/m2)
203	! zcondicea(ngrid,nlayermx): condensation rate in layer l (kg/m2/s)
204	! zcondices(ngrid) : condensation rate on the ground (kg/m2/s)
205	! zfallice(ngrid,nlayermx) : amount of ice falling from layer l (kg/m2/s)
206	! zdtlatent(ngrid,nlayermx): dT/dt due to phase changes (K/s)
207
208	! Tendencies initially set to 0
209	zcondices(1:ngrid) = 0.
210	pdtsrfc(1:ngrid) = 0.
211	pdpsrf(1:ngrid) = 0.
212	ztsrfhist(1:ngrid) = 0.
213	condsub(1:ngrid) = .false.
214	pdicen2(1:ngrid) = 0.
215	zfallheat=0
216	pdqc(1:ngrid,1:nlayer,1:nq)=0
217	pdtc(1:ngrid,1:nlayer)=0
218	pduc(1:ngrid,1:nlayer)=0
219	pdvc(1:ngrid,1:nlayer)=0
220	zfallice(1:ngrid,1:nlayer+1)=0
221	zcondicea(1:ngrid,1:nlayer)=0
222	zdtlatent(1:ngrid,1:nlayer)=0
223	zt(1:ngrid,1:nlayer)=0.
224
225	!-----------------------------------------------------------------------
226	! Atmospheric condensation
227
228	! Condensation / sublimation in the atmosphere
229	! --------------------------------------------
230	! (calcul of zcondicea , zfallice and pdtc)
231
232	zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)+ pdt(1:ngrid,1:nlayer)*ptimestep
233	if (igcm_n2.ne.0) then
234	zqn2(1:ngrid,1:nlayer) = 1. ! & temporaire
235	! zqn2(1:ngrid,1:nlayer)= pq(1:ngrid,1:nlayer,igcm_n2) + pdq(1:ngrid,1:nlayer,igcm_n2)*ptimestep
236	else
237	zqn2(1:ngrid,1:nlayer) = 1.
238	end if
239
240	if (.not.fast) then
241	! Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2
242	DO l=1,nlayer
243	DO ig=1,ngrid
244	ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l))
245	ENDDO
246	ENDDO
247
248	DO l=nlayer,1,-1
249	DO ig=1,ngrid
250	pdtc(ig,l)=0. ! final tendancy temperature set to 0
251
252	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
253	condsub(ig)=.true. !condensation at level l
254	IF (zfallice(ig,l+1).gt.0) then
255	zfallheat=zfallice(ig,l+1)*&
256	(pphi(ig,l+1)-pphi(ig,l) +&
257	cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2
258	ELSE
259	zfallheat=0.
260	ENDIF
261	zdtlatent(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
262	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))&
263	ccondzdtlatent(ig,l)- zfallheat
264	! Case when the ice from above sublimes entirely
265	! """""""""""""""""""""""""""""""""""""""""""""""
266	IF ((zfallice(ig,l+1).lt.-zcondicea(ig,l)) &
267	.AND. (zfallice(ig,l+1).gt.0)) THEN
268
269	zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/&
270	(ccond*(pplev(ig,l)-pplev(ig,l+1)))
271	zcondicea(ig,l)= -zfallice(ig,l+1)
272	END IF
273
274	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
275
276	END IF
277
278	ENDDO
279	ENDDO
280	endif ! not fast
281
282	!-----------------------------------------------------------------------
283	! Condensation/sublimation on the ground
284	! (calculation of zcondices and pdtsrfc)
285
286	! Adding subtimesteps : in fast version, pressures too low lead to
287	! instabilities.
288	IF (fast) THEN
289	IF (pplev(1,1).gt.0.3) THEN
290	nsubtimestep= 1
291	ELSE
292	nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1)
293	ENDIF
294	ELSE
295	nsubtimestep= 1 ! if more then kp and p00 have to be calculated
296	! for nofast mode
297	ENDIF
298	subtimestep=ptimestep/float(nsubtimestep)
299
300	do itsub=1,nsubtimestep
301	! first loop : getting zplev, ztsurf
302	IF (itsub.eq.1) then
303	DO ig=1,ngrid
304	zplev(ig)=pplev(ig,1)
305	ztsrfhist(ig)=ptsrf(ig) + pdtsrf(ig)*ptimestep
306	ztsrf(ig)=ptsrf(ig) + pdtsrf(ig)*subtimestep !!
307	zpicen2(ig)=picen2(ig)
308	ENDDO
309	ELSE
310	! additional loop :
311	! 1) pressure updated
312	! 2) direct redistribution of pressure over the globe
313	! 3) modification pressure for unstable cases
314	! 4) pressure update to conserve mass
315	! 5) temperature updated with radiative tendancies
316	DO ig=1,ngrid
317	zplevhist(ig)=zplev(ig)
318	zplevnew(ig)=zplev(ig)+pdpsrf(ig)*subtimestep ! 1)
319	!IF (zplevnew(ig).le.0.0001) then
320	! zplevnew(ig)=0.0001*kp(ig)/p00
321	!ENDIF
322	ENDDO
323	! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code)
324	globzplevnew=glob_average2d(zplevnew)
325	DO ig=1,ngrid
326	zplev(ig)=kp(ig)*globzplevnew/p00 ! 2)
327	ENDDO
328	DO ig=1,ngrid ! 3) unstable case condensation and sublimation: zplev=zplevhist
329	IF (((pdpsrf(ig).lt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).le.ztsrf(ig))).or. &
330	((pdpsrf(ig).gt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).ge.ztsrf(ig)))) then
331	zplev(ig)=zplevhist(ig)
332	ENDIF
333	zplevhist(ig)=zplev(ig)
334	ENDDO
335	zplev=zplev*globzplevnew/glob_average2d(zplevhist) ! 4)
336	DO ig=1,ngrid ! 5)
337	zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))(ptsrf(ig)4-ztsrf(ig)*4)
338	ztsrf(ig)=ztsrf(ig)+pdtsrfc(ig)subtimestep+zdtsrf(ig)subtimestep
339	zpicen2(ig)=zpicen2(ig)+pdicen2(ig)*subtimestep
340	ENDDO
341	ENDIF ! (itsub=1)
342
343	DO ig=1,ngrid
344	! forecast of frost temperature ztcondsol
345	!ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
346	ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
347
348	! Loop over where we have condensation / sublimation
349	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground cond
350	((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublim
351	(zpicen2(ig) .GT. 0.))) THEN
352	condsub(ig) = .true. ! condensation or sublimation
353
354	! Condensation or partial sublimation of N2 ice
355	if (ztsrf(ig) .LT. ztcondsol(ig)) then ! condensation
356	! Include a correction to account for the cooling of air near
357	! the surface before condensing:
358	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
359	/((lw_n2+cpp(zt(ig,1)-ztcondsol(ig)))subtimestep)
360	else ! sublimation
361	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
362	/(lw_n2*subtimestep)
363	end if
364
365	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep
366
367	! partial sublimation of N2 ice
368	! If the entire N_2 ice layer sublimes
369	! (including what has just condensed in the atmosphere)
370	IF((zpicen2(ig)/subtimestep).LE. &
371	-zcondices(ig))THEN
372	zcondices(ig) = -zpicen2(ig)/subtimestep
373	pdtsrfc(ig)=(lw_n2/pcapcal(ig))* &
374	(zcondices(ig))
375	END IF
376
377	! Changing N2 ice amount and pressure
378
379	pdicen2(ig) = zcondices(ig)
380	pdpsrf(ig) = -pdicen2(ig)*g
381	! pdpsrf(ig) = 0. ! OPTION to check impact N2 sub/cond
382	IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then
383	pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
384	pdicen2(ig)=-pdpsrf(ig)/g
385	ENDIF
386
387	ELSE ! no condsub
388	pdpsrf(ig)=0.
389	pdicen2(ig)=0.
390	pdtsrfc(ig)=0.
391	ENDIF
392	ENDDO ! ig
393	enddo ! subtimestep
394
395	! Updating pressure, temperature and ice reservoir
396	DO ig=1,ngrid
397	pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep
398	! Two options here : 1 ok, 2 is wrong
399	pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep
400	!pdicen2(ig)=-pdpsrf(ig)/g
401
402	pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep
403
404	! security
405	if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then
406	write(,) 'WARNING in condense_n2:'
407	write(,) picen2(ig),pdicen2(ig)*ptimestep
408	pdicen2(ig)= -picen2(ig)/ptimestep
409	pdpsrf(ig)=-pdicen2(ig)*g
410	endif
411
412	if(.not.picen2(ig).ge.0.) THEN
413	! if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
414	print, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)ptimestep
415	! pdicen2(ig)= -picen2(ig)/ptimestep
416	! else
417	picen2(ig)=0.0
418	! endif
419	endif
420	ENDDO
421
422	! ***************************************************************
423	! Correction to account for redistribution between sigma or hybrid
424	! layers when changing surface pressure (and warming/cooling
425	! of the n2 currently changing phase).
426	! *************************************************************
427	if (.not.fast) then
428	DO ig=1,ngrid
429	if (condsub(ig)) then
430
431	! Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
432	! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
433	zmflux(1) = -zcondices(ig)
434	DO l=1,nlayer
435	zmflux(l+1) = zmflux(l) -zcondicea(ig,l) &
436	+ (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
437	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
438	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
439	END DO
440
441	! Mass of each layer
442	! ------------------
443	DO l=1,nlayer
444	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
445	END DO
446
447
448	! Corresponding fluxes for T,U,V,Q
449	! """"""""""""""""""""""""""""""""
450	! averaging operator for TRANSPORT
451	! """"""""""""""""""""""""""""""""
452
453	! Subtimestep loop to perform the redistribution gently and simultaneously with
454	! the other tendencies
455	! Estimation of subtimestep (using only the first layer, the most critical)
456	dtmax=ptimestep
457	if (zmflux(1).gt.1.e-20) then
458	dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
459	endif
460	nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
461	subtimestep=ptimestep/float(nsubtimestep)
462
463	! New flux for each subtimestep
464	do l=1,nlayer+1
465	w(l)=-zmflux(l)*subtimestep
466	enddo
467	! initializing variables that will vary during subtimestep:
468	do l=1,nlayer
469	ztc(l) =pt(ig,l)
470	zu(l) =pu(ig,l)
471	zv(l) =pv(ig,l)
472	do iq=1,nqmx
473	zq(l,iq) = pq(ig,l,iq)
474	enddo
475	end do
476
477	! loop over nsubtimestep
478	! """"""""""""""""""""""
479	do itsub=1,nsubtimestep
480	! Progressively adding tendancies from other processes.
481	do l=1,nlayer
482	ztc(l) =ztc(l) +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
483	zu(l) =zu(l) +pdu( ig,l) * subtimestep
484	zv(l) =zv(l) +pdv( ig,l) * subtimestep
485	do iq=1,nqmx
486	zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
487	enddo
488	end do
489
490	! Change of mass in each layer
491	do l=1,nlayer
492	masse(l)=masse(l)+pdpsrf(ig)subtimestep(pplev(ig,l) - pplev(ig,l+1))&
493	/(g*pplev(ig,1))
494	end do
495
496	ztm(2:nlayermx+1)=0.
497	zum(2:nlayermx+1)=0.
498	zvm(2:nlayermx+1)=0.
499	zqm1(1:nlayermx+1)=0.
500
501	! Van Leer scheme:
502	call vl1d(nlayer,ztc,2.,masse,w,ztm)
503	call vl1d(nlayer,zu ,2.,masse,w,zum)
504	call vl1d(nlayer,zv ,2.,masse,w,zvm)
505	do iq=1,nqmx
506	do l=1,nlayer
507	zq1(l)=zq(l,iq)
508	enddo
509	zqm1(1)=zqm(1,iq)
510	call vl1d(nlayer,zq1,2.,masse,w,zqm1)
511	do l=2,nlayer
512	zqm(l,iq)=zqm1(l)
513	enddo
514	enddo
515
516	! Correction to prevent negative value for qn2
517	if (igcm_n2.ne.0) then
518	zqm(1,igcm_n2)=1.
519	do l=1,nlayer-1
520	if (w(l)zqm(l,igcm_n2).gt.zq(l,igcm_n2)masse(l)) then
521	zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2), &
522	(zqm(l,igcm_n2)w(l) -zq(l,igcm_n2)masse(l))/w(l+1) )
523	else
524	exit
525	endif
526	end do
527	end if
528
529	! Value transfert at the surface interface when condensation sublimation:
530
531	if (zmflux(1).lt.0) then
532	! Surface condensation
533	zum(1)= zu(1)
534	zvm(1)= zv(1)
535	ztm(1) = ztc(1)
536	else
537	! Surface sublimation:
538	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
539	zum(1) = 0
540	zvm(1) = 0
541	end if
542	do iq=1,nqmx
543	zqm(1,iq)=0. ! most tracer do not condense !
544	enddo
545	! Special case if the tracer is n2 gas
546	if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
547
548	!!! Source haze: 0.02 pourcent when n2 sublimes
549	IF (source_haze) THEN
550	IF (pdicen2(ig).lt.0) THEN
551	DO iq=1,nq
552	tname=noms(iq)
553	if (tname(1:4).eq."haze") then
554	!zqm(1,iq)=0.02
555	!zqm(1,iq)=-pdicen2(ig)*0.02
556	zqm(1,iq)=-pdicen2(ig)ptimestep0.02
557	!zqm(10,iq)=-pdicen2(ig)ptimestep100.
558	!zqm(1,iq)=-pdicen2(ig)*1000000.
559
560	endif
561	ENDDO
562	ENDIF
563	ENDIF
564	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
565	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
566	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
567	do iq=1,nqmx
568	zqm(nlayer+1,iq)= zq(nlayer,iq)
569	enddo
570
571	! Tendencies on T, U, V, Q
572	! """""""""""""""""""""""
573	DO l=1,nlayer
574
575	! Tendencies on T
576	zdtsig(ig,l) = (1/masse(l)) * &
577	( zmflux(l)*(ztm(l) - ztc(l)) &
578	- zmflux(l+1)*(ztm(l+1) - ztc(l)) &
579	+ zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
580
581	! Tendencies on U
582	pduc(ig,l) = (1/masse(l)) * &
583	( zmflux(l)*(zum(l) - zu(l))&
584	- zmflux(l+1)*(zum(l+1) - zu(l)) )
585
586	! Tendencies on V
587	pdvc(ig,l) = (1/masse(l)) * &
588	( zmflux(l)*(zvm(l) - zv(l)) &
589	- zmflux(l+1)*(zvm(l+1) - zv(l)) )
590
591	END DO
592
593	! Tendencies on Q
594	do iq=1,nqmx
595	if (iq.eq.igcm_n2) then
596	! SPECIAL Case when the tracer IS N2 :
597	DO l=1,nlayer
598	pdqc(ig,l,iq)= (1/masse(l)) * &
599	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
600	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
601	+ zcondicea(ig,l)*(zq(l,iq)-1.) )
602	END DO
603	else
604	DO l=1,nlayer
605	pdqc(ig,l,iq)= (1/masse(l)) * &
606	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
607	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) &
608	+ zcondicea(ig,l)*zq(l,iq) )
609	END DO
610	end if
611	enddo
612	! Update variables at the end of each subtimestep.
613	do l=1,nlayer
614	ztc(l) =ztc(l) + zdtsig(ig,l) *subtimestep
615	zu(l) =zu(l) + pduc(ig,l) *subtimestep
616	zv(l) =zv(l) + pdvc(ig,l) *subtimestep
617	do iq=1,nqmx
618	zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
619	enddo
620	end do
621	enddo ! loop on nsubtimestep
622	! Recomputing Total tendencies
623	do l=1,nlayer
624	pdtc(ig,l) = (ztc(l) - zt(ig,l) )/ptimestep
625	pduc(ig,l) = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
626	pdvc(ig,l) = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
627	do iq=1,nqmx
628	pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep
629
630
631	! Correction temporaire
632	if (iq.eq.igcm_n2) then
633	if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
634	.lt.0.01) then ! if n2 < 1 % !
635	pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
636	end if
637	end if
638
639	enddo
640	end do
641	! *****************************TEMPORAIRE ****************
642	if (ngrid.eq.1) then
643	write(,) 'nsubtimestep=' ,nsubtimestep
644	write(,) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
645	write(,) 'masse apres' , masse(1)
646	write(,) 'zmfluxDT, l=1' , zmflux(1)ptimestep
647	write(,) 'zmfluxDT, l=2' , zmflux(2)ptimestep
648	write(,) 'pq, l=1,2,3' , pq(1,1,1), pq(1,2,1),pq(1,3,1)
649	write(,) 'zq, l=1,2,3' , zq(1,1), zq(2,1),zq(3,1)
650	write(,) 'dqDt l=1' , pdq(1,1,1)ptimestep
651	write(,) 'dqcondDt l=1' , pdqc(1,1,1)ptimestep
652	end if
653
654	! ***********************************************************
655	end if ! if (condsub)
656	END DO ! loop on ig
657	endif ! not fast
658
659	! ********** Option Olkin to prevent N2 effect in the south******
660	112 continue
661	if (olkin) then
662	DO ig=1,ngrid
663	if (lati(ig).lt.0) then
664	pdtsrfc(ig) = max(0.,pdtsrfc(ig))
665	pdpsrf(ig) = 0.
666	pdicen2(ig) = 0.
667	do l=1,nlayer
668	pdtc(ig,l) = max(0.,zdtlatent(ig,l))
669	pduc(ig,l) = 0.
670	pdvc(ig,l) = 0.
671	do iq=1,nqmx
672	pdqc(ig,l,iq) = 0
673	enddo
674	end do
675	end if
676	END DO
677	end if
678	! *******************************************************************
679
680	! ***************************************************************
681	! Ecriture des diagnostiques
682	! ***************************************************************
683
684	! DO l=1,nlayer
685	! DO ig=1,ngrid
686	! Taux de cond en kg.m-2.pa-1.s-1
687	! tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
688	! Taux de cond en kg.m-3.s-1
689	! tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
690	! END DO
691	! END DO
692	! call WRITEDIAGFI(ngridmx,'tconda1', &
693	! 'Taux de condensation N2 atmospherique /Pa', &
694	! 'kg.m-2.Pa-1.s-1',3,tconda1)
695	! call WRITEDIAGFI(ngridmx,'tconda2', &
696	! 'Taux de condensation N2 atmospherique /m', &
697	! 'kg.m-3.s-1',3,tconda2)
698
699
700	return
701	end subroutine condens_n2
702
703	!-------------------------------------------------------------------------
704
705	real function tcond_n2(p,vmr)
706	! Calculates the condensation temperature for N2 at pressure P and vmr
707	implicit none
708	real, intent(in):: p,vmr
709
710	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
711	IF (p.ge.0.529995) then
712	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
713	! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e50.98)log(1./(0.1255701.e5)p*vmr))
714	tcond_n2 = (1.)/(0.01583606505-1.211938776e-3log(7.963685594e-5p*vmr))
715	ELSE
716	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
717	! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e51.09)log(1./(0.508059)pvmr))
718	tcond_n2 = (1.)/(1./35.6-1.089633028e-3log(1.968275338p*vmr))
719	ENDIF
720	return
721	end function tcond_n2
722
723	!-------------------------------------------------------------------------
724
725	real function pcond_n2(t,vmr)
726	! Calculates the condensation pressure for N2 at temperature T and vmr
727	implicit none
728	real, intent(in):: t,vmr
729
730	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
731	IF (t.ge.35.6) then
732	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
733	! pcond_n2 = 0.1255701.e5/vmrexp((2.5e50.98)/296.925(1./63.1470-1./t))
734	pcond_n2 = 0.125570e5/vmrexp(825.1241896(1./63.147-1./t))
735	ELSE
736	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB
737	! pcond_n2 = 0.508059/vmrexp((2.5e51.09)/296.925*(1./35.6-1./t))
738	pcond_n2 = 0.508059/vmrexp(917.7401701(1./35.6-1./t))
739	ENDIF
740	return
741	end function pcond_n2
742
743	!-------------------------------------------------------------------------
744
745	real function glob_average2d(var)
746	! Calculates the global average of variable var
747	use comgeomfi_h
748	implicit none
749	#include "dimensions.h"
750	#include "dimphys.h"
751
752	! INTEGER ngrid
753	REAL var(ngridmx)
754	INTEGER ig
755
756	glob_average2d = 0.
757	DO ig=2,ngridmx-1
758	glob_average2d = glob_average2d + var(ig)*area(ig)
759	END DO
760	glob_average2d = glob_average2d + var(1)area(1)iim
761	glob_average2d = glob_average2d + var(ngridmx)area(ngridmx)iim
762	glob_average2d = glob_average2d/(totarea+(area(1)+area(ngridmx))*(iim-1))
763
764	end function glob_average2d
765
766	! *****************************************************************
767
768	subroutine vl1d(nlayer,q,pente_max,masse,w,qm)
769	!
770	! Operateur de moyenne inter-couche pour calcul de transport type
771	! Van-Leer " pseudo amont " dans la verticale
772	! q,w sont des arguments d'entree pour le s-pg ....
773	! masse : masse de la couche Dp/g
774	! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
775	! pente_max = 2 conseillee
776	! --------------------------------------------------------------------
777	IMPLICIT NONE
778
779	#include "dimensions.h"
780
781	! Arguments:
782	! ----------
783	integer nlayer
784	real masse(nlayer),pente_max
785	REAL q(nlayer),qm(nlayer+1)
786	REAL w(nlayer+1)
787	!
788	! Local
789	! ---------
790	!
791	INTEGER l
792	!
793	real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
794	real sigw, Mtot, MQtot
795	integer m
796
797
798	! On oriente tout dans le sens de la pression
799	! W > 0 WHEN DOWN !!!!!!!!!!!!!
800
801	do l=2,nlayer
802	dzqw(l)=q(l-1)-q(l)
803	adzqw(l)=abs(dzqw(l))
804	enddo
805
806	do l=2,nlayer-1
807	if(dzqw(l)*dzqw(l+1).gt.0.) then
808	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
809	else
810	dzq(l)=0.
811	endif
812	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
813	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
814	enddo
815
816	dzq(1)=0.
817	dzq(nlayer)=0.
818
819	do l = 1,nlayer-1
820
821	! Regular scheme (transfered mass < layer mass)
822	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
823	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
824	sigw=w(l+1)/masse(l+1)
825	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
826	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
827	sigw=w(l+1)/masse(l)
828	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
829
830	! Extended scheme (transfered mass > layer mass)
831	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
832	else if(w(l+1).gt.0.) then
833	m=l+1
834	Mtot = masse(m)
835	MQtot = masse(m)*q(m)
836	!if (m.lt.nlayer) then ! because some compilers will have problems
837	! ! evaluating masse(nlayer+1)
838	do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
839	m=m+1
840	Mtot = Mtot + masse(m)
841	MQtot = MQtot + masse(m)*q(m)
842	! if (m.eq.nlayer) exit
843	end do
844	!endif
845	if (m.lt.nlayer) then
846	sigw=(w(l+1)-Mtot)/masse(m+1)
847	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot) &
848	(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
849	else
850	w(l+1) = Mtot
851	qm(l+1) = Mqtot / Mtot
852	write(,) 'top layer is disapearing !'
853	stop
854	end if
855	else ! if(w(l+1).lt.0)
856	m = l-1
857	Mtot = masse(m+1)
858	MQtot = masse(m+1)*q(m+1)
859	if (m.gt.0) then ! because some compilers will have problems
860	! evaluating masse(0)
861	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
862	m=m-1
863	Mtot = Mtot + masse(m+1)
864	MQtot = MQtot + masse(m+1)*q(m+1)
865	if (m.eq.0) exit
866	end do
867	endif
868	if (m.gt.0) then
869	sigw=(w(l+1)+Mtot)/masse(m)
870	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot) &
871	(q(m)-0.5(1.+sigw)dzq(m)) )
872	else
873	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
874	end if
875	end if
876	enddo
877
878	return
879	end subroutine vl1d
880

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

Download in other formats:

Original Format