Context Navigation

condense_n2.F90 @ 3917

Last change on this file since 3917 was 3910, checked in by tbertrand, 4 months ago

Pluto PCM:

adding option to use N2 ice fractional maps (n2frac) read in the startfi.nc
adding option in newstart to correct (tsurf and tsoil) for too warm or too cold N2-free (correct_t_non2) or N2-rich (correct_t_n2) patches

Property svn:executable set to *

File size: 27.8 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,n2frac)
7
8	use comgeomfi_h
9	use comcstfi_mod, only: g, r, cpp, pi
10	USE surfdat_h, only: phisfi,kp,p00
11	USE tracer_h, only: noms, igcm_n2, lw_n2
12	USE callkeys_mod, only: fast,ch4lag,latlag,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag
13	USE vertical_layers_mod, ONLY: ap,bp
14	USE geometry_mod, only: latitude, longitude, cell_area
15	use planetwide_mod, only: planetwide_sumval
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	REAL n2frac(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,SAVE :: firstcall=.true.
149	!$OMP THREADPRIVATE(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	zcondices(ig) = zcondices(ig)*n2frac(ig) !N2frac added as ponderation
375
376	! partial sublimation of N2 ice
377	! If the entire N_2 ice layer sublimes
378	! (including what has just condensed in the atmosphere)
379	IF((zpicen2(ig)/subtimestep).LE. &
380	-zcondices(ig))THEN
381	zcondices(ig) = -zpicen2(ig)/subtimestep
382	pdtsrfc(ig)=(lw_n2/pcapcal(ig))* &
383	(zcondices(ig))
384	END IF
385
386	! Changing N2 ice amount and pressure
387	pdicen2(ig) = zcondices(ig)
388	pdpsrf(ig) = -pdicen2(ig)*g
389	! pdpsrf(ig) = 0. ! OPTION to check impact N2 sub/cond
390	IF (fast.and.(zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001)) then
391	pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
392	pdicen2(ig)=-pdpsrf(ig)/g
393	ENDIF
394
395	ELSE ! no condsub
396
397	!IF (zpicen2(ig) .GT. 0.) THEN
398	! print,'TB24 no condsub',latitude(ig)180./pi,longitude(ig)*180./pi
399	!ENDIF
400	pdpsrf(ig)=0.
401	pdicen2(ig)=0.
402	pdtsrfc(ig)=0.
403	ENDIF
404	ENDDO ! ig
405	enddo ! itsub,nsubtimestep
406
407	! Updating pressure, temperature and ice reservoir
408	DO ig=1,klon
409	pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep
410	! Two options here : 1 ok, 2 is wrong
411	pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep
412	!pdicen2(ig)=-pdpsrf(ig)/g
413
414	pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep
415
416	!IF (picen2(ig) .GT. 0. .and. abs(pdicen2(ig)).lt.1.e-12) THEN
417	! print,'TB24 low pdq',latitude(ig)180./pi,longitude(ig)*180./pi
418	!ENDIF
419
420
421
422	! security
423	if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then
424	write(,) 'WARNING in condense_n2:'
425	write(,) picen2(ig),pdicen2(ig)*ptimestep
426	pdicen2(ig)= -picen2(ig)/ptimestep
427	pdpsrf(ig)=-pdicen2(ig)*g
428	endif
429
430	if(.not.picen2(ig).ge.0.) THEN
431	! if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
432	print, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)ptimestep
433	! pdicen2(ig)= -picen2(ig)/ptimestep
434	! else
435	picen2(ig)=0.0
436	! endif
437	endif
438	ENDDO ! ig
439
440	! ***************************************************************
441	! Correction to account for redistribution between sigma or hybrid
442	! layers when changing surface pressure (and warming/cooling
443	! of the n2 currently changing phase).
444	! *************************************************************
445	if (.not.fast) then
446	DO ig=1,klon
447	if (condsub(ig)) then
448
449	! Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
450	! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
451	zmflux(1) = -zcondices(ig)
452	DO l=1,klev
453	zmflux(l+1) = zmflux(l) -zcondicea(ig,l) &
454	+ (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
455	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
456	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
457	END DO
458
459	! Mass of each layer
460	! ------------------
461	DO l=1,klev
462	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
463	END DO
464
465	! Corresponding fluxes for T,U,V,Q
466	! """"""""""""""""""""""""""""""""
467	! averaging operator for TRANSPORT
468	! """"""""""""""""""""""""""""""""
469
470	! Subtimestep loop to perform the redistribution gently and simultaneously with
471	! the other tendencies
472	! Estimation of subtimestep (using only the first layer, the most critical)
473	dtmax=ptimestep
474	if (zmflux(1).gt.1.e-20) then
475	dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
476	endif
477	nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
478	subtimestep=ptimestep/float(nsubtimestep)
479
480	! New flux for each subtimestep
481	do l=1,klev+1
482	w(l)=-zmflux(l)*subtimestep
483	enddo
484	! initializing variables that will vary during subtimestep:
485	do l=1,klev
486	ztc(l) =pt(ig,l)
487	zu(l) =pu(ig,l)
488	zv(l) =pv(ig,l)
489	do iq=1,nq
490	zq(l,iq) = pq(ig,l,iq)
491	enddo
492	end do
493
494	! loop over nsubtimestep
495	! """"""""""""""""""""""
496	do itsub=1,nsubtimestep
497	! Progressively adding tendancies from other processes.
498	do l=1,klev
499	ztc(l) =ztc(l) +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
500	zu(l) =zu(l) +pdu( ig,l) * subtimestep
501	zv(l) =zv(l) +pdv( ig,l) * subtimestep
502	do iq=1,nq
503	zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
504	enddo
505	end do
506
507	! Change of mass in each layer
508	do l=1,klev
509	masse(l)=masse(l)+pdpsrf(ig)subtimestep(pplev(ig,l) - pplev(ig,l+1))&
510	/(g*pplev(ig,1))
511	end do
512
513	ztm(:)=0.
514	zum(:)=0.
515	zvm(:)=0.
516	zqm1(:)=0.
517	zqm(:,:)=0.
518
519	! Van Leer scheme:
520	call vl1d(klev,ztc,2.,masse,w,ztm)
521	call vl1d(klev,zu ,2.,masse,w,zum)
522	call vl1d(klev,zv ,2.,masse,w,zvm)
523	do iq=1,nq
524	do l=1,klev
525	zq1(l)=zq(l,iq)
526	enddo
527	zqm1(1)=zqm(1,iq)
528	call vl1d(klev,zq1,2.,masse,w,zqm1)
529	do l=2,klev
530	zqm(l,iq)=zqm1(l)
531	enddo
532	enddo
533
534	! Correction to prevent negative value for qn2
535	if (igcm_n2.ne.0) then
536	zqm(1,igcm_n2)=1.
537	do l=1,klev-1
538	if (w(l)zqm(l,igcm_n2).gt.zq(l,igcm_n2)masse(l)) then
539	zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2), &
540	(zqm(l,igcm_n2)w(l) -zq(l,igcm_n2)masse(l))/w(l+1) )
541	else
542	exit
543	endif
544	end do
545	end if
546
547	! Value transfert at the surface interface when condensation sublimation:
548	if (zmflux(1).lt.0) then
549	! Surface condensation
550	zum(1)= zu(1)
551	zvm(1)= zv(1)
552	ztm(1) = ztc(1)
553	else
554	! Surface sublimation:
555	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
556	zum(1) = 0
557	zvm(1) = 0
558	end if
559	do iq=1,nq
560	zqm(1,iq)=0. ! most tracer do not condense !
561	enddo
562	! Special case if the tracer is n2 gas
563	if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
564
565	!!! Source haze: 0.02 pourcent when n2 sublimes
566	IF (source_haze) THEN
567	IF (pdicen2(ig).lt.0) THEN
568	DO iq=1,nq
569	tname=noms(iq)
570	if (tname(1:4).eq."haze") then
571	!zqm(1,iq)=0.02
572	!zqm(1,iq)=-pdicen2(ig)*0.02
573	zqm(1,iq)=-pdicen2(ig)ptimestep0.02
574	!zqm(10,iq)=-pdicen2(ig)ptimestep100.
575	!zqm(1,iq)=-pdicen2(ig)*1000000.
576
577	endif
578	ENDDO
579	ENDIF
580	ENDIF
581	ztm(klev+1)= ztc(klev) ! should not be used, but...
582	zum(klev+1)= zu(klev) ! should not be used, but...
583	zvm(klev+1)= zv(klev) ! should not be used, but...
584	do iq=1,nq
585	zqm(klev+1,iq)= zq(klev,iq)
586	enddo
587
588	! Tendencies on T, U, V, Q
589	! """""""""""""""""""""""
590	DO l=1,klev
591
592	! Tendencies on T
593	zdtsig(ig,l) = (1/masse(l)) * &
594	( zmflux(l)*(ztm(l) - ztc(l)) &
595	- zmflux(l+1)*(ztm(l+1) - ztc(l)) &
596	+ zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
597
598	! Tendencies on U
599	pduc(ig,l) = (1/masse(l)) * &
600	( zmflux(l)*(zum(l) - zu(l))&
601	- zmflux(l+1)*(zum(l+1) - zu(l)) )
602
603	! Tendencies on V
604	pdvc(ig,l) = (1/masse(l)) * &
605	( zmflux(l)*(zvm(l) - zv(l)) &
606	- zmflux(l+1)*(zvm(l+1) - zv(l)) )
607
608	END DO
609
610	! Tendencies on Q
611	do iq=1,nq
612	if (iq.eq.igcm_n2) then
613	! SPECIAL Case when the tracer IS N2 :
614	DO l=1,klev
615	pdqc(ig,l,iq)= (1/masse(l)) * &
616	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
617	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
618	+ zcondicea(ig,l)*(zq(l,iq)-1.) )
619	END DO
620	else
621	DO l=1,klev
622	pdqc(ig,l,iq)= (1/masse(l)) * &
623	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
624	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) &
625	+ zcondicea(ig,l)*zq(l,iq) )
626	END DO
627	end if
628	enddo
629	! Update variables at the end of each subtimestep.
630	do l=1,klev
631	ztc(l) =ztc(l) + zdtsig(ig,l) *subtimestep
632	zu(l) =zu(l) + pduc(ig,l) *subtimestep
633	zv(l) =zv(l) + pdvc(ig,l) *subtimestep
634	do iq=1,nq
635	zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
636	enddo
637	end do
638
639	enddo ! loop on nsubtimestep
640
641	! Recomputing Total tendencies
642	do l=1,klev
643	pdtc(ig,l) = (ztc(l) - zt(ig,l) )/ptimestep
644	pduc(ig,l) = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
645	pdvc(ig,l) = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
646	do iq=1,nq
647	pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep
648
649
650	! Correction temporaire
651	if (iq.eq.igcm_n2) then
652	if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
653	.lt.0.01) then ! if n2 < 1 % !
654	write(,) 'Warning: n2 < 1%'
655	pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
656	end if
657	end if
658
659	enddo
660	end do
661	! *****************************TEMPORAIRE ****************
662	if (klon.eq.1) then
663	write(,) 'nsubtimestep=' ,nsubtimestep
664	write(,) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
665	write(,) 'masse apres' , masse(1)
666	write(,) 'zmfluxDT, l=1' , zmflux(1)ptimestep
667	write(,) 'zmfluxDT, l=2' , zmflux(2)ptimestep
668	write(,) 'pq, l=1,2,3' , pq(1,1,1), pq(1,2,1),pq(1,3,1)
669	write(,) 'zq, l=1,2,3' , zq(1,1), zq(2,1),zq(3,1)
670	write(,) 'dqDt l=1' , pdq(1,1,1)ptimestep
671	write(,) 'dqcondDt l=1' , pdqc(1,1,1)ptimestep
672	end if
673
674	! ***********************************************************
675	end if ! if (condsub)
676	END DO ! loop on ig
677	endif ! not fast
678
679	! ***************************************************************
680	! Ecriture des diagnostiques
681	! ***************************************************************
682
683	! DO l=1,klev
684	! DO ig=1,klon
685	! Taux de cond en kg.m-2.pa-1.s-1
686	! tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
687	! Taux de cond en kg.m-3.s-1
688	! tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
689	! END DO
690	! END DO
691	! call WRITEDIAGFI(klon,'tconda1', &
692	! 'Taux de condensation N2 atmospherique /Pa', &
693	! 'kg.m-2.Pa-1.s-1',3,tconda1)
694	! call WRITEDIAGFI(klon,'tconda2', &
695	! 'Taux de condensation N2 atmospherique /m', &
696	! 'kg.m-3.s-1',3,tconda2)
697
698
699	return
700	end subroutine condense_n2
701
702	!-------------------------------------------------------------------------
703
704	real function tcond_n2(p,vmr)
705	! Calculates the condensation temperature for N2 at pressure P and vmr
706	implicit none
707	real, intent(in):: p,vmr
708
709	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
710	IF (p.ge.0.529995) then
711	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
712	! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e50.98)log(1./(0.1255701.e5)p*vmr))
713	tcond_n2 = (1.)/(0.01583606505-1.211938776e-3log(7.963685594e-5p*vmr))
714	ELSE
715	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
716	! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e51.09)log(1./(0.508059)pvmr))
717	tcond_n2 = (1.)/(1./35.6-1.089633028e-3log(1.968275338p*vmr))
718	ENDIF
719	return
720	end function tcond_n2
721
722	!-------------------------------------------------------------------------
723
724	real function pcond_n2(t,vmr)
725	! Calculates the condensation pressure for N2 at temperature T and vmr
726	implicit none
727	real, intent(in):: t,vmr
728
729	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
730	IF (t.ge.35.6) then
731	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
732	! pcond_n2 = 0.1255701.e5/vmrexp((2.5e50.98)/296.925(1./63.1470-1./t))
733	pcond_n2 = 0.125570e5/vmrexp(825.1241896(1./63.147-1./t))
734	ELSE
735	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB
736	! pcond_n2 = 0.508059/vmrexp((2.5e51.09)/296.925*(1./35.6-1./t))
737	pcond_n2 = 0.508059/vmrexp(917.7401701(1./35.6-1./t))
738	ENDIF
739	return
740	end function pcond_n2
741
742	!-------------------------------------------------------------------------
743
744	subroutine vl1d(klev,q,pente_max,masse,w,qm)
745	!
746	! Operateur de moyenne inter-couche pour calcul de transport type
747	! Van-Leer " pseudo amont " dans la verticale
748	! q,w sont des arguments d'entree pour le s-pg ....
749	! masse : masse de la couche Dp/g
750	! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
751	! pente_max = 2 conseillee
752	! --------------------------------------------------------------------
753	IMPLICIT NONE
754
755	! Arguments:
756	! ----------
757	integer klev
758	real masse(klev),pente_max
759	REAL q(klev),qm(klev+1)
760	REAL w(klev+1)
761	!
762	! Local
763	! ---------
764	!
765	INTEGER l
766	!
767	real dzq(klev),dzqw(klev),adzqw(klev),dzqmax
768	real sigw, Mtot, MQtot
769	integer m
770
771
772	! On oriente tout dans le sens de la pression
773	! W > 0 WHEN DOWN !!!!!!!!!!!!!
774
775	do l=2,klev
776	dzqw(l)=q(l-1)-q(l)
777	adzqw(l)=abs(dzqw(l))
778	enddo
779
780	do l=2,klev-1
781	if(dzqw(l)*dzqw(l+1).gt.0.) then
782	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
783	else
784	dzq(l)=0.
785	endif
786	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
787	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
788	enddo
789
790	dzq(1)=0.
791	dzq(klev)=0.
792
793	do l = 1,klev-1
794
795	! Regular scheme (transfered mass < layer mass)
796	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
797	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
798	sigw=w(l+1)/masse(l+1)
799	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
800	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
801	sigw=w(l+1)/masse(l)
802	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
803
804	! Extended scheme (transfered mass > layer mass)
805	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
806	else if(w(l+1).gt.0.) then
807	m=l+1
808	Mtot = masse(m)
809	MQtot = masse(m)*q(m)
810	if (m.lt.klev) then ! because some compilers will have problems
811	! evaluating masse(klev+1)
812	do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1))))
813	m=m+1
814	Mtot = Mtot + masse(m)
815	MQtot = MQtot + masse(m)*q(m)
816	if (m.eq.klev) exit
817	end do
818	endif
819	if (m.lt.klev) then
820	sigw=(w(l+1)-Mtot)/masse(m+1)
821	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot) &
822	(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
823	else
824	w(l+1) = Mtot
825	qm(l+1) = Mqtot / Mtot
826	print*,"m", m
827	print*,"Mtot", Mtot
828	print*,"Mqtot", Mqtot
829	call abort_physic("condense_n2","top layer is disappearing !",1)
830	end if
831	else ! if(w(l+1).lt.0)
832	m = l-1
833	Mtot = masse(m+1)
834	MQtot = masse(m+1)*q(m+1)
835	if (m.gt.0) then ! because some compilers will have problems
836	! evaluating masse(0)
837	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
838	m=m-1
839	Mtot = Mtot + masse(m+1)
840	MQtot = MQtot + masse(m+1)*q(m+1)
841	if (m.eq.0) exit
842	end do
843	endif
844	if (m.gt.0) then
845	sigw=(w(l+1)+Mtot)/masse(m)
846	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot) &
847	(q(m)-0.5(1.+sigw)dzq(m)) )
848	else
849	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
850	end if
851	end if
852	enddo
853
854	return
855	end subroutine vl1d

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

Download in other formats:

Original Format