Context Navigation

condense_n2.F90 @ 3558

Last change on this file since 3558 was 3539, checked in by tbertrand, 5 weeks ago
LMDZ.PLUTO: Update and Validation of the Volatile Transport Model (VTM, or "nogcm") Merging missing elements with Pluto.old TB
Property svn:executable set to ``*
File size: 27.6 KB

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

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

Download in other formats:

Original Format