Context Navigation

condense_n2.F90 @ 3380

Last change on this file since 3380 was 3232, checked in by afalco, 9 months ago
Pluto PCM: 1D functional AF
Property svn:executable set to ``*
File size: 28.3 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 comvert_mod, ONLY: ap,bp
15	use geometry_mod, only: latitude
16
17
18	implicit none
19
20	!==================================================================
21	! Purpose
22	! -------
23	! Condense and/or sublime N2 ice on the ground and in the
24	! atmosphere, and sediment the ice.
25	!
26	! Inputs
27	! ------
28	! klon Number of vertical columns
29	! klev Number of layers
30	! pplay(klon,klev) Pressure layers
31	! pplev(klon,klev+1) Pressure levels
32	! pt(klon,klev) Temperature (in K)
33	! ptsrf(klon) Surface temperature
34	!
35	! pdt(klon,klev) Time derivative before condensation/sublimation of pt
36	! pdtsrf(klon) Time derivative before condensation/sublimation of ptsrf
37	! picen2(klon) n2 ice at the surface (kg/m2)
38	!
39	! Outputs
40	! -------
41	! pdpsrf(klon) \ Contribution of condensation/sublimation
42	! pdtc(klon,klev) / to the time derivatives of Ps, pt, and ptsrf
43	! pdtsrfc(klon) /
44	! pdicen2(klon) Tendancy n2 ice at the surface (kg/m2)
45	!
46	! Both
47	! ----
48	! psolaralb(klon) Albedo at the surface
49	! pemisurf(klon) Emissivity of the surface
50	!
51	! Authors
52	! -------
53	! Francois Forget (1996,2013)
54	! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
55	!
56	!
57	!==================================================================
58
59	!-----------------------------------------------------------------------
60	! Arguments
61
62	INTEGER klon, klev, nq
63
64	REAL ptimestep
65	REAL pplay(klon,klev),pplev(klon,klev+1)
66	REAL pcapcal(klon)
67	REAL pt(klon,klev)
68	REAL ptsrf(klon),flu1(klon),flu2(klon),flu3(klon)
69	REAL pphi(klon,klev)
70	REAL pdt(klon,klev),pdtsrf(klon),pdtc(klon,klev)
71	REAL pdtsrfc(klon),pdpsrf(klon)
72	REAL picen2(klon),psolaralb(klon),pemisurf(klon)
73
74
75	REAL pu(klon,klev) , pv(klon,klev)
76	REAL pdu(klon,klev) , pdv(klon,klev)
77	REAL pduc(klon,klev) , pdvc(klon,klev)
78	REAL pq(klon,klev,nq),pdq(klon,klev,nq)
79	REAL pdqc(klon,klev,nq)
80
81	!-----------------------------------------------------------------------
82	! Local variables
83
84	INTEGER l,ig,ilay,it,iq,ich4_gas
85
86	REAL*8 zt(klon,klev)
87	REAL tcond_n2
88	REAL pcond_n2
89	REAL glob_average2d ! function
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 zplevnew(klon)
126	REAL zplev(klon)
127	REAL zpicen2(klon)
128	REAL ztsrfhist(klon)
129	REAL zdtsrf(klon)
130	REAL globzplevnew
131
132	! REAL vmrn2(klon)
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 condense_n2cloud: ccond=',ccond,' latcond=',lw_n2
171
172	! calculate global mean surface pressure for the fast mode
173	IF (fast) THEN
174	IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon))
175	DO ig=1,klon
176	kp(ig) = exp(-phisfi(ig)/(r*38.))
177	ENDDO
178	p00=glob_average2d(kp) ! mean pres at ref level
179	ENDIF
180
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
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	globzplevnew=glob_average2d(zplevnew)
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	zplev=zplev*globzplevnew/glob_average2d(zplevhist) ! 4)
337	DO ig=1,klon ! 5)
338	zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))(ptsrf(ig)4-ztsrf(ig)*4)
339	ztsrf(ig)=ztsrf(ig)+pdtsrfc(ig)subtimestep+zdtsrf(ig)subtimestep
340	zpicen2(ig)=zpicen2(ig)+pdicen2(ig)*subtimestep
341	ENDDO
342	ENDIF ! (itsub=1)
343
344	DO ig=1,klon
345	! forecast of frost temperature ztcondsol
346	ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
347	!ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
348
349	! Loop over where we have condensation / sublimation
350	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground cond
351	((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublim
352	(zpicen2(ig) .GT. 0.))) THEN
353	condsub(ig) = .true. ! condensation or sublimation
354
355	! Condensation or partial sublimation of N2 ice
356	if (ztsrf(ig) .LT. ztcondsol(ig)) then ! condensation
357	! Include a correction to account for the cooling of air near
358	! the surface before condensing:
359	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
360	/((lw_n2+cpp(zt(ig,1)-ztcondsol(ig)))subtimestep)
361	else ! sublimation
362	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
363	/(lw_n2*subtimestep)
364	end if
365
366	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep
367
368	! partial sublimation of N2 ice
369	! If the entire N_2 ice layer sublimes
370	! (including what has just condensed in the atmosphere)
371	IF((zpicen2(ig)/subtimestep).LE. &
372	-zcondices(ig))THEN
373	zcondices(ig) = -zpicen2(ig)/subtimestep
374	pdtsrfc(ig)=(lw_n2/pcapcal(ig))* &
375	(zcondices(ig))
376	END IF
377
378	! Changing N2 ice amount and pressure
379
380	pdicen2(ig) = zcondices(ig)
381	pdpsrf(ig) = -pdicen2(ig)*g
382	! pdpsrf(ig) = 0. ! OPTION to check impact N2 sub/cond
383	IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then
384	pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
385	pdicen2(ig)=-pdpsrf(ig)/g
386	ENDIF
387
388	ELSE ! no condsub
389	pdpsrf(ig)=0.
390	pdicen2(ig)=0.
391	pdtsrfc(ig)=0.
392	ENDIF
393	ENDDO ! ig
394	enddo ! subtimestep
395
396	! Updating pressure, temperature and ice reservoir
397	DO ig=1,klon
398	pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep
399	! Two options here : 1 ok, 2 is wrong
400	pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep
401	!pdicen2(ig)=-pdpsrf(ig)/g
402
403	pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep
404
405	! security
406	if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then
407	write(,) 'WARNING in condense_n2:'
408	write(,) picen2(ig),pdicen2(ig)*ptimestep
409	pdicen2(ig)= -picen2(ig)/ptimestep
410	pdpsrf(ig)=-pdicen2(ig)*g
411	endif
412
413	if(.not.picen2(ig).ge.0.) THEN
414	! if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
415	print, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)ptimestep
416	! pdicen2(ig)= -picen2(ig)/ptimestep
417	! else
418	picen2(ig)=0.0
419	! endif
420	endif
421	ENDDO
422
423	! ***************************************************************
424	! Correction to account for redistribution between sigma or hybrid
425	! layers when changing surface pressure (and warming/cooling
426	! of the n2 currently changing phase).
427	! *************************************************************
428	if (.not.fast) then
429	DO ig=1,klon
430	if (condsub(ig)) then
431
432	! Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
433	! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
434	zmflux(1) = -zcondices(ig)
435	DO l=1,klev
436	zmflux(l+1) = zmflux(l) -zcondicea(ig,l) &
437	+ (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
438	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
439	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
440	END DO
441
442	! Mass of each layer
443	! ------------------
444	DO l=1,klev
445	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
446	END DO
447
448
449	! Corresponding fluxes for T,U,V,Q
450	! """"""""""""""""""""""""""""""""
451	! averaging operator for TRANSPORT
452	! """"""""""""""""""""""""""""""""
453
454	! Subtimestep loop to perform the redistribution gently and simultaneously with
455	! the other tendencies
456	! Estimation of subtimestep (using only the first layer, the most critical)
457	dtmax=ptimestep
458	if (zmflux(1).gt.1.e-20) then
459	dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
460	endif
461	nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
462	subtimestep=ptimestep/float(nsubtimestep)
463
464	! New flux for each subtimestep
465	do l=1,klev+1
466	w(l)=-zmflux(l)*subtimestep
467	enddo
468	! initializing variables that will vary during subtimestep:
469	do l=1,klev
470	ztc(l) =pt(ig,l)
471	zu(l) =pu(ig,l)
472	zv(l) =pv(ig,l)
473	do iq=1,nq
474	zq(l,iq) = pq(ig,l,iq)
475	enddo
476	end do
477
478	! loop over nsubtimestep
479	! """"""""""""""""""""""
480	do itsub=1,nsubtimestep
481	! Progressively adding tendancies from other processes.
482	do l=1,klev
483	ztc(l) =ztc(l) +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
484	zu(l) =zu(l) +pdu( ig,l) * subtimestep
485	zv(l) =zv(l) +pdv( ig,l) * subtimestep
486	do iq=1,nq
487	zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
488	enddo
489	end do
490
491	! Change of mass in each layer
492	do l=1,klev
493	masse(l)=masse(l)+pdpsrf(ig)subtimestep(pplev(ig,l) - pplev(ig,l+1))&
494	/(g*pplev(ig,1))
495	end do
496
497	! Value transfert at the surface interface when condensation sublimation:
498
499	if (zmflux(1).lt.0) then
500	! Surface condensation
501	zum(1)= zu(1)
502	zvm(1)= zv(1)
503	ztm(1) = ztc(1)
504	else
505	! Surface sublimation:
506	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
507	zum(1) = 0
508	zvm(1) = 0
509	end if
510	do iq=1,nq
511	zqm(1,iq)=0. ! most tracer do not condense !
512	enddo
513	! Special case if the tracer is n2 gas
514	if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
515
516	ztm(2:klev+1)=0.
517	zum(2:klev+1)=0.
518	zvm(2:klev+1)=0.
519	zqm1(1:klev+1)=0.
520
521	! Van Leer scheme:
522	call vl1d(klev,ztc,2.,masse,w,ztm)
523	call vl1d(klev,zu ,2.,masse,w,zum)
524	call vl1d(klev,zv ,2.,masse,w,zvm)
525	do iq=1,nq
526	do l=1,klev
527	zq1(l)=zq(l,iq)
528	enddo
529	zqm1(1)=zqm(1,iq)
530	call vl1d(klev,zq1,2.,masse,w,zqm1)
531	do l=2,klev
532	zqm(l,iq)=zqm1(l)
533	enddo
534	enddo
535
536	! Correction to prevent negative value for qn2
537	if (igcm_n2.ne.0) then
538	zqm(1,igcm_n2)=1.
539	do l=1,klev-1
540	if (w(l)zqm(l,igcm_n2).gt.zq(l,igcm_n2)masse(l)) then
541	zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2), &
542	(zqm(l,igcm_n2)w(l) -zq(l,igcm_n2)masse(l))/w(l+1) )
543	else
544	exit
545	endif
546	end do
547	end if
548
549
550	!!! Source haze: 0.02 pourcent when n2 sublimes
551	IF (source_haze) THEN
552	IF (pdicen2(ig).lt.0) THEN
553	DO iq=1,nq
554	tname=noms(iq)
555	if (tname(1:4).eq."haze") then
556	!zqm(1,iq)=0.02
557	!zqm(1,iq)=-pdicen2(ig)*0.02
558	zqm(1,iq)=-pdicen2(ig)ptimestep0.02
559	!zqm(10,iq)=-pdicen2(ig)ptimestep100.
560	!zqm(1,iq)=-pdicen2(ig)*1000000.
561
562	endif
563	ENDDO
564	ENDIF
565	ENDIF
566	ztm(klev+1)= ztc(klev) ! should not be used, but...
567	zum(klev+1)= zu(klev) ! should not be used, but...
568	zvm(klev+1)= zv(klev) ! should not be used, but...
569	do iq=1,nq
570	zqm(klev+1,iq)= zq(klev,iq)
571	enddo
572
573	! Tendencies on T, U, V, Q
574	! """""""""""""""""""""""
575	DO l=1,klev
576
577	! Tendencies on T
578	zdtsig(ig,l) = (1/masse(l)) * &
579	( zmflux(l)*(ztm(l) - ztc(l)) &
580	- zmflux(l+1)*(ztm(l+1) - ztc(l)) &
581	+ zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
582
583	! Tendencies on U
584	pduc(ig,l) = (1/masse(l)) * &
585	( zmflux(l)*(zum(l) - zu(l))&
586	- zmflux(l+1)*(zum(l+1) - zu(l)) )
587
588	! Tendencies on V
589	pdvc(ig,l) = (1/masse(l)) * &
590	( zmflux(l)*(zvm(l) - zv(l)) &
591	- zmflux(l+1)*(zvm(l+1) - zv(l)) )
592
593	END DO
594
595	! Tendencies on Q
596	do iq=1,nq
597	if (iq.eq.igcm_n2) then
598	! SPECIAL Case when the tracer IS N2 :
599	DO l=1,klev
600	pdqc(ig,l,iq)= (1/masse(l)) * &
601	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
602	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
603	+ zcondicea(ig,l)*(zq(l,iq)-1.) )
604	END DO
605	else
606	DO l=1,klev
607	pdqc(ig,l,iq)= (1/masse(l)) * &
608	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
609	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) &
610	+ zcondicea(ig,l)*zq(l,iq) )
611	END DO
612	end if
613	enddo
614	! Update variables at the end of each subtimestep.
615	do l=1,klev
616	ztc(l) =ztc(l) + zdtsig(ig,l) *subtimestep
617	zu(l) =zu(l) + pduc(ig,l) *subtimestep
618	zv(l) =zv(l) + pdvc(ig,l) *subtimestep
619	do iq=1,nq
620	zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
621	enddo
622	end do
623	enddo ! loop on nsubtimestep
624	! Recomputing Total tendencies
625	do l=1,klev
626	pdtc(ig,l) = (ztc(l) - zt(ig,l) )/ptimestep
627	pduc(ig,l) = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
628	pdvc(ig,l) = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
629	do iq=1,nq
630	pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep
631
632
633	! Correction temporaire
634	if (iq.eq.igcm_n2) then
635	if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
636	.lt.0.01) then ! if n2 < 1 % !
637	write(,) 'Warning: n2 < 1%'
638	pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
639	end if
640	end if
641
642	enddo
643	end do
644	! *****************************TEMPORAIRE ****************
645	if (klon.eq.1) then
646	write(,) 'nsubtimestep=' ,nsubtimestep
647	write(,) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
648	write(,) 'masse apres' , masse(1)
649	write(,) 'zmfluxDT, l=1' , zmflux(1)ptimestep
650	write(,) 'zmfluxDT, l=2' , zmflux(2)ptimestep
651	write(,) 'pq, l=1,2,3' , pq(1,1,1), pq(1,2,1),pq(1,3,1)
652	write(,) 'zq, l=1,2,3' , zq(1,1), zq(2,1),zq(3,1)
653	write(,) 'dqDt l=1' , pdq(1,1,1)ptimestep
654	write(,) 'dqcondDt l=1' , pdqc(1,1,1)ptimestep
655	end if
656
657	! ***********************************************************
658	end if ! if (condsub)
659	END DO ! loop on ig
660	endif ! not fast
661
662	! ********** Option Olkin to prevent N2 effect in the south******
663	112 continue
664	if (olkin) then
665	DO ig=1,klon
666	if (latitude(ig).lt.0) then
667	pdtsrfc(ig) = max(0.,pdtsrfc(ig))
668	pdpsrf(ig) = 0.
669	pdicen2(ig) = 0.
670	do l=1,klev
671	pdtc(ig,l) = max(0.,zdtlatent(ig,l))
672	pduc(ig,l) = 0.
673	pdvc(ig,l) = 0.
674	do iq=1,nq
675	pdqc(ig,l,iq) = 0
676	enddo
677	end do
678	end if
679	END DO
680	end if
681	! *******************************************************************
682
683	! ***************************************************************
684	! Ecriture des diagnostiques
685	! ***************************************************************
686
687	! DO l=1,klev
688	! DO ig=1,klon
689	! Taux de cond en kg.m-2.pa-1.s-1
690	! tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
691	! Taux de cond en kg.m-3.s-1
692	! tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
693	! END DO
694	! END DO
695	! call WRITEDIAGFI(klon,'tconda1', &
696	! 'Taux de condensation N2 atmospherique /Pa', &
697	! 'kg.m-2.Pa-1.s-1',3,tconda1)
698	! call WRITEDIAGFI(klon,'tconda2', &
699	! 'Taux de condensation N2 atmospherique /m', &
700	! 'kg.m-3.s-1',3,tconda2)
701
702
703	return
704	end subroutine condense_n2
705
706	!-------------------------------------------------------------------------
707
708	real function tcond_n2(p,vmr)
709	! Calculates the condensation temperature for N2 at pressure P and vmr
710	implicit none
711	real, intent(in):: p,vmr
712
713	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
714	IF (p.ge.0.529995) then
715	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
716	! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e50.98)log(1./(0.1255701.e5)p*vmr))
717	tcond_n2 = (1.)/(0.01583606505-1.211938776e-3log(7.963685594e-5p*vmr))
718	ELSE
719	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
720	! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e51.09)log(1./(0.508059)pvmr))
721	tcond_n2 = (1.)/(1./35.6-1.089633028e-3log(1.968275338p*vmr))
722	ENDIF
723	return
724	end function tcond_n2
725
726	!-------------------------------------------------------------------------
727
728	real function pcond_n2(t,vmr)
729	! Calculates the condensation pressure for N2 at temperature T and vmr
730	implicit none
731	real, intent(in):: t,vmr
732
733	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
734	IF (t.ge.35.6) then
735	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
736	! pcond_n2 = 0.1255701.e5/vmrexp((2.5e50.98)/296.925(1./63.1470-1./t))
737	pcond_n2 = 0.125570e5/vmrexp(825.1241896(1./63.147-1./t))
738	ELSE
739	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB
740	! pcond_n2 = 0.508059/vmrexp((2.5e51.09)/296.925*(1./35.6-1./t))
741	pcond_n2 = 0.508059/vmrexp(917.7401701(1./35.6-1./t))
742	ENDIF
743	return
744	end function pcond_n2
745
746	!-------------------------------------------------------------------------
747
748	real function glob_average2d(var)
749	! Calculates the global average of variable var
750	use comgeomfi_h
751	use dimphy, only: klon
752	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
753	use geometry_mod, only: cell_area, latitude
754
755	implicit none
756
757	! INTEGER klon
758	REAL var(klon)
759	INTEGER ig
760
761	glob_average2d = 0.
762	DO ig=2,klon-1
763	glob_average2d = glob_average2d + var(ig)*cell_area(ig)
764	END DO
765	glob_average2d = glob_average2d + var(1)cell_area(1)nbp_lon
766	glob_average2d = glob_average2d + var(klon)cell_area(klon)nbp_lon
767	glob_average2d = glob_average2d/(totarea+(cell_area(1)+cell_area(klon))*(nbp_lon-1))
768
769	end function glob_average2d
770
771	! *****************************************************************
772
773	subroutine vl1d(klev,q,pente_max,masse,w,qm)
774	!
775	! Operateur de moyenne inter-couche pour calcul de transport type
776	! Van-Leer " pseudo amont " dans la verticale
777	! q,w sont des arguments d'entree pour le s-pg ....
778	! masse : masse de la couche Dp/g
779	! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
780	! pente_max = 2 conseillee
781	! --------------------------------------------------------------------
782	IMPLICIT NONE
783
784	! Arguments:
785	! ----------
786	integer klev
787	real masse(klev),pente_max
788	REAL q(klev),qm(klev+1)
789	REAL w(klev+1)
790	!
791	! Local
792	! ---------
793	!
794	INTEGER l
795	!
796	real dzq(klev),dzqw(klev),adzqw(klev),dzqmax
797	real sigw, Mtot, MQtot
798	integer m
799
800
801	! On oriente tout dans le sens de la pression
802	! W > 0 WHEN DOWN !!!!!!!!!!!!!
803
804	do l=2,klev
805	dzqw(l)=q(l-1)-q(l)
806	adzqw(l)=abs(dzqw(l))
807	enddo
808
809	do l=2,klev-1
810	if(dzqw(l)*dzqw(l+1).gt.0.) then
811	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
812	else
813	dzq(l)=0.
814	endif
815	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
816	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
817	enddo
818
819	dzq(1)=0.
820	dzq(klev)=0.
821
822	do l = 1,klev-1
823
824	! Regular scheme (transfered mass < layer mass)
825	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
826	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
827	sigw=w(l+1)/masse(l+1)
828	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
829	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
830	sigw=w(l+1)/masse(l)
831	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
832
833	! Extended scheme (transfered mass > layer mass)
834	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
835	else if(w(l+1).gt.0.) then
836	m=l+1
837	Mtot = masse(m)
838	MQtot = masse(m)*q(m)
839	!if (m.lt.klev) then ! because some compilers will have problems
840	! ! evaluating masse(klev+1)
841	do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1))))
842	m=m+1
843	Mtot = Mtot + masse(m)
844	MQtot = MQtot + masse(m)*q(m)
845	! if (m.eq.klev) exit
846	end do
847	!endif
848	if (m.lt.klev) then
849	sigw=(w(l+1)-Mtot)/masse(m+1)
850	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot) &
851	(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
852	else
853	w(l+1) = Mtot
854	qm(l+1) = Mqtot / Mtot
855	write(,) 'top layer is disapearing !'
856	stop
857	end if
858	else ! if(w(l+1).lt.0)
859	m = l-1
860	Mtot = masse(m+1)
861	MQtot = masse(m+1)*q(m+1)
862	if (m.gt.0) then ! because some compilers will have problems
863	! evaluating masse(0)
864	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
865	m=m-1
866	Mtot = Mtot + masse(m+1)
867	MQtot = MQtot + masse(m+1)*q(m+1)
868	if (m.eq.0) exit
869	end do
870	endif
871	if (m.gt.0) then
872	sigw=(w(l+1)+Mtot)/masse(m)
873	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot) &
874	(q(m)-0.5(1.+sigw)dzq(m)) )
875	else
876	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
877	end if
878	end if
879	enddo
880
881	return
882	end subroutine vl1d
883

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

Download in other formats:

Original Format