Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

condens_n2.F90_saved @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 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.2 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	firstcall=.false.
189	ENDIF
190
191	!-----------------------------------------------------------------------
192	! Calculation of n2 condensation / sublimation
193
194	! Variables used:
195	! picen2(ngrid) : amount of n2 ice on the ground (kg/m2)
196	! zcondicea(ngrid,nlayermx): condensation rate in layer l (kg/m2/s)
197	! zcondices(ngrid) : condensation rate on the ground (kg/m2/s)
198	! zfallice(ngrid,nlayermx) : amount of ice falling from layer l (kg/m2/s)
199	! zdtlatent(ngrid,nlayermx): dT/dt due to phase changes (K/s)
200
201	! Tendencies initially set to 0
202	zcondices(1:ngrid) = 0.
203	pdtsrfc(1:ngrid) = 0.
204	pdpsrf(1:ngrid) = 0.
205	ztsrfhist(1:ngrid) = 0.
206	condsub(1:ngrid) = .false.
207	pdicen2(1:ngrid) = 0.
208	zfallheat=0
209	pdqc(1:ngrid,1:nlayer,1:nq)=0
210	pdtc(1:ngrid,1:nlayer)=0
211	pduc(1:ngrid,1:nlayer)=0
212	pdvc(1:ngrid,1:nlayer)=0
213	zfallice(1:ngrid,1:nlayer+1)=0
214	zcondicea(1:ngrid,1:nlayer)=0
215	zdtlatent(1:ngrid,1:nlayer)=0
216	zt(1:ngrid,1:nlayer)=0.
217
218	!-----------------------------------------------------------------------
219	! Atmospheric condensation
220
221	! Condensation / sublimation in the atmosphere
222	! --------------------------------------------
223	! (calcul of zcondicea , zfallice and pdtc)
224
225	zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)+ pdt(1:ngrid,1:nlayer)*ptimestep
226	if (igcm_n2.ne.0) then
227	zqn2(1:ngrid,1:nlayer) = 1. ! & temporaire
228	! zqn2(1:ngrid,1:nlayer)= pq(1:ngrid,1:nlayer,igcm_n2) + pdq(1:ngrid,1:nlayer,igcm_n2)*ptimestep
229	else
230	zqn2(1:ngrid,1:nlayer) = 1.
231	end if
232
233	if (.not.fast) then
234	! Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2
235	DO l=1,nlayer
236	DO ig=1,ngrid
237	ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l))
238	ENDDO
239	ENDDO
240
241	DO l=nlayer,1,-1
242	DO ig=1,ngrid
243	pdtc(ig,l)=0. ! final tendancy temperature set to 0
244
245	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
246	condsub(ig)=.true. !condensation at level l
247	IF (zfallice(ig,l+1).gt.0) then
248	zfallheat=zfallice(ig,l+1)*&
249	(pphi(ig,l+1)-pphi(ig,l) +&
250	cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2
251	ELSE
252	zfallheat=0.
253	ENDIF
254	zdtlatent(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
255	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))&
256	ccondzdtlatent(ig,l)- zfallheat
257	! Case when the ice from above sublimes entirely
258	! """""""""""""""""""""""""""""""""""""""""""""""
259	IF ((zfallice(ig,l+1).lt.-zcondicea(ig,l)) &
260	.AND. (zfallice(ig,l+1).gt.0)) THEN
261
262	zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/&
263	(ccond*(pplev(ig,l)-pplev(ig,l+1)))
264	zcondicea(ig,l)= -zfallice(ig,l+1)
265	END IF
266
267	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
268
269	END IF
270
271	ENDDO
272	ENDDO
273	endif ! not fast
274
275	!-----------------------------------------------------------------------
276	! Condensation/sublimation on the ground
277	! (calculation of zcondices and pdtsrfc)
278
279	! Adding subtimesteps : in fast version, pressures too low lead to
280	! instabilities.
281	IF (fast) THEN
282	IF (pplev(1,1).gt.0.3) THEN
283	nsubtimestep= 1
284	ELSE
285	nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1)
286	ENDIF
287	ELSE
288	nsubtimestep= 1 ! if more then kp and p00 have to be calculated
289	! for nofast mode
290	ENDIF
291	subtimestep=ptimestep/float(nsubtimestep)
292
293	do itsub=1,nsubtimestep
294	! first loop : getting zplev, ztsurf
295	IF (itsub.eq.1) then
296	DO ig=1,ngrid
297	zplev(ig)=pplev(ig,1)
298	ztsrfhist(ig)=ptsrf(ig) + pdtsrf(ig)*ptimestep
299	ztsrf(ig)=ptsrf(ig) + pdtsrf(ig)*subtimestep !!
300	zpicen2(ig)=picen2(ig)
301	ENDDO
302	ELSE
303	! additional loop :
304	! 1) pressure updated
305	! 2) direct redistribution of pressure over the globe
306	! 3) modification pressure for unstable cases
307	! 4) pressure update to conserve mass
308	! 5) temperature updated with radiative tendancies
309	DO ig=1,ngrid
310	zplevhist(ig)=zplev(ig)
311	zplevnew(ig)=zplev(ig)+pdpsrf(ig)*subtimestep ! 1)
312	!IF (zplevnew(ig).le.0.0001) then
313	! zplevnew(ig)=0.0001*kp(ig)/p00
314	!ENDIF
315	ENDDO
316	! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code)
317	globzplevnew=glob_average2d(zplevnew)
318	DO ig=1,ngrid
319	zplev(ig)=kp(ig)*globzplevnew/p00 ! 2)
320	ENDDO
321	DO ig=1,ngrid ! 3) unstable case condensation and sublimation: zplev=zplevhist
322	IF (((pdpsrf(ig).lt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).le.ztsrf(ig))).or. &
323	((pdpsrf(ig).gt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).ge.ztsrf(ig)))) then
324	zplev(ig)=zplevhist(ig)
325	ENDIF
326	zplevhist(ig)=zplev(ig)
327	ENDDO
328	zplev=zplev*globzplevnew/glob_average2d(zplevhist) ! 4)
329	DO ig=1,ngrid ! 5)
330	zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))(ptsrf(ig)4-ztsrf(ig)*4)
331	ztsrf(ig)=ztsrf(ig)+pdtsrfc(ig)subtimestep+zdtsrf(ig)subtimestep
332	zpicen2(ig)=zpicen2(ig)+pdicen2(ig)*subtimestep
333	ENDDO
334	ENDIF ! (itsub=1)
335
336	DO ig=1,ngrid
337	! forecast of frost temperature ztcondsol
338	!ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
339	ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
340
341	! Loop over where we have condensation / sublimation
342	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground cond
343	((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublim
344	(zpicen2(ig) .GT. 0.))) THEN
345	condsub(ig) = .true. ! condensation or sublimation
346
347	! Condensation or partial sublimation of N2 ice
348	if (ztsrf(ig) .LT. ztcondsol(ig)) then ! condensation
349	! Include a correction to account for the cooling of air near
350	! the surface before condensing:
351	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
352	/((lw_n2+cpp(zt(ig,1)-ztcondsol(ig)))subtimestep)
353	else ! sublimation
354	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
355	/(lw_n2*subtimestep)
356	end if
357
358	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep
359
360	! partial sublimation of N2 ice
361	! If the entire N_2 ice layer sublimes
362	! (including what has just condensed in the atmosphere)
363	IF((zpicen2(ig)/subtimestep).LE. &
364	-zcondices(ig))THEN
365	zcondices(ig) = -zpicen2(ig)/subtimestep
366	pdtsrfc(ig)=(lw_n2/pcapcal(ig))* &
367	(zcondices(ig))
368	END IF
369
370	! Changing N2 ice amount and pressure
371
372	pdicen2(ig) = zcondices(ig)
373	pdpsrf(ig) = -pdicen2(ig)*g
374	! pdpsrf(ig) = 0. ! OPTION to check impact N2 sub/cond
375	IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then
376	pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
377	pdicen2(ig)=-pdpsrf(ig)/g
378	ENDIF
379
380	ELSE ! no condsub
381	pdpsrf(ig)=0.
382	pdicen2(ig)=0.
383	pdtsrfc(ig)=0.
384	ENDIF
385	ENDDO ! ig
386	enddo ! subtimestep
387
388	! Updating pressure, temperature and ice reservoir
389	DO ig=1,ngrid
390	pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep
391	! Two options here : 1 ok, 2 is wrong
392	pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep
393	!pdicen2(ig)=-pdpsrf(ig)/g
394
395	pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep
396
397	! security
398	if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then
399	write(,) 'WARNING in condense_n2:'
400	write(,) picen2(ig),pdicen2(ig)*ptimestep
401	pdicen2(ig)= -picen2(ig)/ptimestep
402	pdpsrf(ig)=-pdicen2(ig)*g
403	endif
404
405	if(.not.picen2(ig).ge.0.) THEN
406	! if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
407	print, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)ptimestep
408	! pdicen2(ig)= -picen2(ig)/ptimestep
409	! else
410	picen2(ig)=0.0
411	! endif
412	endif
413	ENDDO
414
415	! ***************************************************************
416	! Correction to account for redistribution between sigma or hybrid
417	! layers when changing surface pressure (and warming/cooling
418	! of the n2 currently changing phase).
419	! *************************************************************
420	if (.not.fast) then
421	DO ig=1,ngrid
422	if (condsub(ig)) then
423
424	! Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
425	! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
426
427	zmflux(1) = -zcondices(ig)
428	DO l=1,nlayer
429	zmflux(l+1) = zmflux(l) -zcondicea(ig,l) &
430	+ (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
431	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
432	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
433	END DO
434
435	! Mass of each layer
436	! ------------------
437	DO l=1,nlayer
438	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
439	END DO
440
441
442	! Corresponding fluxes for T,U,V,Q
443	! """"""""""""""""""""""""""""""""
444	! averaging operator for TRANSPORT
445	! """"""""""""""""""""""""""""""""
446
447	! Subtimestep loop to perform the redistribution gently and simultaneously with
448	! the other tendencies
449	! Estimation of subtimestep (using only the first layer, the most critical)
450	dtmax=ptimestep
451	if (zmflux(1).gt.1.e-20) then
452	dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
453	endif
454	nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
455	subtimestep=ptimestep/float(nsubtimestep)
456
457	! New flux for each subtimestep
458	do l=1,nlayer+1
459	w(l)=-zmflux(l)*subtimestep
460	enddo
461	! initializing variables that will vary during subtimestep:
462	do l=1,nlayer
463	ztc(l) =pt(ig,l)
464	zu(l) =pu(ig,l)
465	zv(l) =pv(ig,l)
466	do iq=1,nqmx
467	zq(l,iq) = pq(ig,l,iq)
468	enddo
469	end do
470
471	! loop over nsubtimestep
472	! """"""""""""""""""""""
473	do itsub=1,nsubtimestep
474	! Progressively adding tendancies from other processes.
475	do l=1,nlayer
476	ztc(l) =ztc(l) +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
477	zu(l) =zu(l) +pdu( ig,l) * subtimestep
478	zv(l) =zv(l) +pdv( ig,l) * subtimestep
479	do iq=1,nqmx
480	zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
481	enddo
482	end do
483
484	! Change of mass in each layer
485	do l=1,nlayer
486	masse(l)=masse(l)+pdpsrf(ig)subtimestep(pplev(ig,l) - pplev(ig,l+1))&
487	/(g*pplev(ig,1))
488	end do
489
490	ztm(2:nlayermx+1)=0.
491	zum(2:nlayermx+1)=0.
492	zvm(2:nlayermx+1)=0.
493	zqm1(1:nlayermx+1)=0.
494
495	! Van Leer scheme:
496	call vl1d(ztc,2.,masse,w,ztm)
497	call vl1d(zu ,2.,masse,w,zum)
498	call vl1d(zv ,2.,masse,w,zvm)
499	do iq=1,nqmx
500	do l=1,nlayer
501	zq1(l)=zq(l,iq)
502	enddo
503	zqm1(1)=zqm(1,iq)
504	call vl1d(zq1,2.,masse,w,zqm1)
505	do l=2,nlayer
506	zqm(l,iq)=zqm1(l)
507	enddo
508	enddo
509
510	! Correction to prevent negative value for qn2
511	if (igcm_n2.ne.0) then
512	do l=1,nlayer-1
513	if (w(l)zqm(l,igcm_n2).gt.zq(l,igcm_n2)masse(l)) then
514	zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2), &
515	(zqm(l,igcm_n2)w(l) -zq(l,igcm_n2)masse(l))/w(l+1) )
516	else
517	exit
518	endif
519	end do
520	end if
521
522	! Value transfert at the surface interface when condensation sublimation:
523
524	if (zmflux(1).lt.0) then
525	! Surface condensation
526	zum(1)= zu(1)
527	zvm(1)= zv(1)
528	ztm(1) = ztc(1)
529	else
530	! Surface sublimation:
531	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
532	zum(1) = 0
533	zvm(1) = 0
534	end if
535	do iq=1,nqmx
536	zqm(1,iq)=0. ! most tracer do not condense !
537	enddo
538	! Special case if the tracer is n2 gas
539	if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
540
541	!!! Source haze: 0.02 pourcent when n2 sublimes
542	IF (source_haze) THEN
543	IF (pdicen2(ig).lt.0) THEN
544	DO iq=1,nq
545	tname=noms(iq)
546	if (tname(1:4).eq."haze") then
547	!zqm(1,iq)=0.02
548	!zqm(1,iq)=-pdicen2(ig)*0.02
549	zqm(1,iq)=-pdicen2(ig)ptimestep0.02
550	!zqm(10,iq)=-pdicen2(ig)ptimestep100.
551	!zqm(1,iq)=-pdicen2(ig)*1000000.
552
553	endif
554	ENDDO
555	ENDIF
556	ENDIF
557	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
558	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
559	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
560	do iq=1,nqmx
561	zqm(nlayer+1,iq)= zq(nlayer,iq)
562	enddo
563
564	! Tendencies on T, U, V, Q
565	! """""""""""""""""""""""
566	DO l=1,nlayer
567
568	! Tendencies on T
569	zdtsig(ig,l) = (1/masse(l)) * &
570	( zmflux(l)*(ztm(l) - ztc(l)) &
571	- zmflux(l+1)*(ztm(l+1) - ztc(l)) &
572	+ zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
573
574	! Tendencies on U
575	pduc(ig,l) = (1/masse(l)) * &
576	( zmflux(l)*(zum(l) - zu(l))&
577	- zmflux(l+1)*(zum(l+1) - zu(l)) )
578
579	! Tendencies on V
580	pdvc(ig,l) = (1/masse(l)) * &
581	( zmflux(l)*(zvm(l) - zv(l)) &
582	- zmflux(l+1)*(zvm(l+1) - zv(l)) )
583
584	END DO
585
586	! Tendencies on Q
587	do iq=1,nqmx
588	if (iq.eq.igcm_n2) then
589	! SPECIAL Case when the tracer IS N2 :
590	DO l=1,nlayer
591	pdqc(ig,l,iq)= (1/masse(l)) * &
592	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
593	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
594	+ zcondicea(ig,l)*(zq(l,iq)-1.) )
595	END DO
596	else
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) )
602	END DO
603	end if
604	enddo
605	! Update variables at the end of each subtimestep.
606	do l=1,nlayer
607	ztc(l) =ztc(l) + zdtsig(ig,l) *subtimestep
608	zu(l) =zu(l) + pduc(ig,l) *subtimestep
609	zv(l) =zv(l) + pdvc(ig,l) *subtimestep
610	do iq=1,nqmx
611	zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
612	enddo
613	end do
614	enddo ! loop on nsubtimestep
615	! Recomputing Total tendencies
616	do l=1,nlayer
617	pdtc(ig,l) = (ztc(l) - zt(ig,l) )/ptimestep
618	pduc(ig,l) = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
619	pdvc(ig,l) = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
620	do iq=1,nqmx
621	pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep
622
623
624	! Correction temporaire
625	if (iq.eq.igcm_n2) then
626	if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
627	.lt.0.01) then ! if n2 < 1 % !
628	pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
629	end if
630	end if
631
632	enddo
633	end do
634	! *****************************TEMPORAIRE ****************
635	if (ngrid.eq.1) then
636	write(,) 'nsubtimestep=' ,nsubtimestep
637	write(,) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
638	write(,) 'masse apres' , masse(1)
639	write(,) 'zmfluxDT, l=1' , zmflux(1)ptimestep
640	write(,) 'zmfluxDT, l=2' , zmflux(2)ptimestep
641	write(,) 'pq, l=1,2,3' , pq(1,1,1), pq(1,2,1),pq(1,3,1)
642	write(,) 'zq, l=1,2,3' , zq(1,1), zq(2,1),zq(3,1)
643	write(,) 'dqDt l=1' , pdq(1,1,1)ptimestep
644	write(,) 'dqcondDt l=1' , pdqc(1,1,1)ptimestep
645	end if
646
647	! ***********************************************************
648	end if ! if (condsub)
649	END DO ! loop on ig
650	endif ! not fast
651
652	! ********** Option Olkin to prevent N2 effect in the south******
653	112 continue
654	if (olkin) then
655	DO ig=1,ngrid
656	if (lati(ig).lt.0) then
657	pdtsrfc(ig) = max(0.,pdtsrfc(ig))
658	pdpsrf(ig) = 0.
659	pdicen2(ig) = 0.
660	do l=1,nlayer
661	pdtc(ig,l) = max(0.,zdtlatent(ig,l))
662	pduc(ig,l) = 0.
663	pdvc(ig,l) = 0.
664	do iq=1,nqmx
665	pdqc(ig,l,iq) = 0
666	enddo
667	end do
668	end if
669	END DO
670	end if
671	! *******************************************************************
672
673	! ***************************************************************
674	! Ecriture des diagnostiques
675	! ***************************************************************
676
677	! DO l=1,nlayer
678	! DO ig=1,ngrid
679	! Taux de cond en kg.m-2.pa-1.s-1
680	! tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
681	! Taux de cond en kg.m-3.s-1
682	! tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
683	! END DO
684	! END DO
685	! call WRITEDIAGFI(ngridmx,'tconda1', &
686	! 'Taux de condensation N2 atmospherique /Pa', &
687	! 'kg.m-2.Pa-1.s-1',3,tconda1)
688	! call WRITEDIAGFI(ngridmx,'tconda2', &
689	! 'Taux de condensation N2 atmospherique /m', &
690	! 'kg.m-3.s-1',3,tconda2)
691
692
693	return
694	end subroutine condens_n2
695
696	!-------------------------------------------------------------------------
697
698	real function tcond_n2(p,vmr)
699	! Calculates the condensation temperature for N2 at pressure P and vmr
700	implicit none
701	real, intent(in):: p,vmr
702
703	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
704	IF (p.ge.0.529995) then
705	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
706	! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e50.98)log(1./(0.1255701.e5)p*vmr))
707	tcond_n2 = (1.)/(0.01583606505-1.211938776e-3log(7.963685594e-5p*vmr))
708	ELSE
709	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
710	! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e51.09)log(1./(0.508059)pvmr))
711	tcond_n2 = (1.)/(1./35.6-1.089633028e-3log(1.968275338p*vmr))
712	ENDIF
713	return
714	end function tcond_n2
715
716	!-------------------------------------------------------------------------
717
718	real function pcond_n2(t,vmr)
719	! Calculates the condensation pressure for N2 at temperature T and vmr
720	implicit none
721	real, intent(in):: t,vmr
722
723	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
724	IF (t.ge.35.6) then
725	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
726	! pcond_n2 = 0.1255701.e5/vmrexp((2.5e50.98)/296.925(1./63.1470-1./t))
727	pcond_n2 = 0.125570e5/vmrexp(825.1241896(1./63.147-1./t))
728	ELSE
729	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB
730	! pcond_n2 = 0.508059/vmrexp((2.5e51.09)/296.925*(1./35.6-1./t))
731	pcond_n2 = 0.508059/vmrexp(917.7401701(1./35.6-1./t))
732	ENDIF
733	return
734	end function pcond_n2
735
736	!-------------------------------------------------------------------------
737
738	real function glob_average2d(var)
739	! Calculates the global average of variable var
740	use comgeomfi_h
741	implicit none
742	#include "dimensions.h"
743	#include "dimphys.h"
744
745	! INTEGER ngrid
746	REAL var(ngridmx)
747	INTEGER ig
748
749	glob_average2d = 0.
750	DO ig=2,ngridmx-1
751	glob_average2d = glob_average2d + var(ig)*area(ig)
752	END DO
753	glob_average2d = glob_average2d + var(1)area(1)iim
754	glob_average2d = glob_average2d + var(ngridmx)area(ngridmx)iim
755	glob_average2d = glob_average2d/(totarea+(area(1)+area(ngridmx))*(iim-1))
756
757	end function glob_average2d
758
759	! *****************************************************************
760
761	subroutine vl1d(q,pente_max,masse,w,qm)
762	!
763	! Operateur de moyenne inter-couche pour calcul de transport type
764	! Van-Leer " pseudo amont " dans la verticale
765	! q,w sont des arguments d'entree pour le s-pg ....
766	! masse : masse de la couche Dp/g
767	! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
768	! pente_max = 2 conseillee
769	! --------------------------------------------------------------------
770	IMPLICIT NONE
771
772	#include "dimensions.h"
773
774	! Arguments:
775	! ----------
776	real masse(llm),pente_max
777	REAL q(llm),qm(llm+1)
778	REAL w(llm+1)
779	!
780	! Local
781	! ---------
782	!
783	INTEGER l
784	!
785	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
786	real sigw, Mtot, MQtot
787	integer m
788
789
790	! On oriente tout dans le sens de la pression
791	! W > 0 WHEN DOWN !!!!!!!!!!!!!
792
793	do l=2,llm
794	dzqw(l)=q(l-1)-q(l)
795	adzqw(l)=abs(dzqw(l))
796	enddo
797
798	do l=2,llm-1
799	if(dzqw(l)*dzqw(l+1).gt.0.) then
800	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
801	else
802	dzq(l)=0.
803	endif
804	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
805	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
806	enddo
807
808	dzq(1)=0.
809	dzq(llm)=0.
810
811	do l = 1,llm-1
812
813	! Regular scheme (transfered mass < layer mass)
814	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
815	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
816	sigw=w(l+1)/masse(l+1)
817	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
818	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
819	sigw=w(l+1)/masse(l)
820	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
821
822	! Extended scheme (transfered mass > layer mass)
823	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
824	else if(w(l+1).gt.0.) then
825	m=l+1
826	Mtot = masse(m)
827	MQtot = masse(m)*q(m)
828	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1))))
829	m=m+1
830	Mtot = Mtot + masse(m)
831	MQtot = MQtot + masse(m)*q(m)
832	end do
833	if (m.lt.llm) then
834	sigw=(w(l+1)-Mtot)/masse(m+1)
835	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot) &
836	(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
837	else
838	w(l+1) = Mtot
839	qm(l+1) = Mqtot / Mtot
840	write(,) 'top layer is disapearing !'
841	stop
842	end if
843	else ! if(w(l+1).lt.0)
844	m = l-1
845	Mtot = masse(m+1)
846	MQtot = masse(m+1)*q(m+1)
847	do while((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(max(m,1)))))
848	m=m-1
849	Mtot = Mtot + masse(m+1)
850	MQtot = MQtot + masse(m+1)*q(m+1)
851	end do
852	if (m.gt.0) then
853	sigw=(w(l+1)+Mtot)/masse(m)
854	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot) &
855	(q(m)-0.5(1.+sigw)dzq(m)) )
856	else
857	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
858	end if
859	end if
860	enddo
861
862	return
863	end subroutine vl1d
864

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

Download in other formats:

Original Format