Context Navigation

condense_n2.F90 @ 3754

Last change on this file since 3754 was 3668, checked in by debatzbr, 4 months ago
Minor correction \ B
Property svn:executable set to ``*
File size: 27.7 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 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
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 zqn2(klon,klev) ! N2 MMR used to compute Tcond/zqn2
90	REAL ztcond (klon,klev)
91	REAL ztcondsol(klon),zfallheat
92	REAL pdicen2(klon)
93	REAL zcondicea(klon,klev), zcondices(klon)
94	REAL zfallice(klon,klev+1)
95	REAL zmflux(klev+1)
96	REAL zu(klev),zv(klev)
97	REAL zq(klev,nq),zq1(klev)
98	REAL ztsrf(klon)
99	REAL ztc(klev), ztm(klev+1)
100	REAL zum(klev+1) , zvm(klev+1)
101	REAL zqm(klev+1,nq),zqm1(klev+1)
102	LOGICAL condsub(klon)
103	REAL subptimestep
104	Integer Ntime
105	real masse (klev),w(klev+1)
106	real wq(klon,klev+1)
107	real vstokes,reff
108	real dWtotsn2
109	real condnconsn2(klon)
110	real nconsMAXn2
111	! Special diagnostic variables
112	real tconda1(klon,klev)
113	real tconda2(klon,klev)
114	real zdtsig (klon,klev)
115	real zdtlatent (klon,klev)
116	real zdt (klon,klev)
117	! REAL albediceF(klon)
118	! SAVE albediceF
119	INTEGER nsubtimestep,itsub !number of subtimestep when calling vl1d
120	REAL subtimestep !ptimestep/nsubtimestep
121	REAL dtmax
122
123	REAL zplevhist(klon)
124	REAL zplevhistglob
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,dimension(:),save,allocatable :: vmrn2
133	!$OMP THREADPRIVATE(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,SAVE :: firstcall=.true.
148	!$OMP THREADPRIVATE(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
160	CHARACTER(len=10) :: tname
161
162	!-----------------------------------------------------------------------
163
164	! Initialisation
165	IF (firstcall) THEN
166	ccond=cpp/(g*lw_n2)
167	print*,'In condense_n2cloud: ccond=',ccond,' latcond=',lw_n2
168
169	! calculate global mean surface pressure for the fast mode
170	IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon))
171	DO ig=1,klon
172	kp(ig) = exp(-phisfi(ig)/(r*38.))
173	ENDDO
174	IF (fast) THEN
175	! mean pres at ref level
176	call planetwide_sumval(kp(:)*cell_area(:)/totarea_planet,p00)
177	ENDIF
178
179	ALLOCATE(vmrn2(klon))
180	vmrn2(:) = 1.
181	!IF (ch4lag) then
182	! DO ig=1,klon
183	! if (latitude(ig)*180./pi.ge.latlag) then
184	! vmrn2(ig) = vmrlag
185	! endif
186	! ENDDO
187	!ENDIF
188	IF (no_n2frost) then
189	DO ig=1,klon
190	if (picen2(ig).eq.0.) then
191	vmrn2(ig) = 1.e-15
192	endif
193	ENDDO
194	ENDIF
195	firstcall=.false.
196	ENDIF
197
198	!-----------------------------------------------------------------------
199	! Calculation of n2 condensation / sublimation
200
201	! Variables used:
202	! picen2(klon) : amount of n2 ice on the ground (kg/m2)
203	! zcondicea(klon,klev): condensation rate in layer l (kg/m2/s)
204	! zcondices(klon) : condensation rate on the ground (kg/m2/s)
205	! zfallice(klon,klev) : amount of ice falling from layer l (kg/m2/s)
206	! zdtlatent(klon,klev): dT/dt due to phase changes (K/s)
207
208	! Tendencies initially set to 0
209	zcondices(1:klon) = 0.
210	pdtsrfc(1:klon) = 0.
211	pdpsrf(1:klon) = 0.
212	ztsrfhist(1:klon) = 0.
213	condsub(1:klon) = .false.
214	pdicen2(1:klon) = 0.
215	zfallheat=0.
216	pdqc(1:klon,1:klev,1:nq)=0.
217	pdtc(1:klon,1:klev)=0.
218	pduc(1:klon,1:klev)=0.
219	pdvc(1:klon,1:klev)=0.
220	zfallice(1:klon,1:klev+1)=0.
221	zcondicea(1:klon,1:klev)=0.
222	zdtlatent(1:klon,1:klev)=0.
223	zt(1:klon,1:klev)=0.
224
225	!-----------------------------------------------------------------------
226	! Atmospheric condensation
227
228	! Condensation / sublimation in the atmosphere
229	! --------------------------------------------
230	! (calcul of zcondicea , zfallice and pdtc)
231
232	zt(1:klon,1:klev)=pt(1:klon,1:klev)+ pdt(1:klon,1:klev)*ptimestep
233	if (igcm_n2.ne.0) then
234	zqn2(1:klon,1:klev) = 1. ! & temporaire
235	! zqn2(1:klon,1:klev)= pq(1:klon,1:klev,igcm_n2) + pdq(1:klon,1:klev,igcm_n2)*ptimestep
236	else
237	zqn2(1:klon,1:klev) = 1.
238	end if
239
240	if (.not.fast) then
241	! Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2
242	DO l=1,klev
243	DO ig=1,klon
244	ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l))
245	ENDDO
246	ENDDO
247
248	DO l=klev,1,-1
249	DO ig=1,klon
250	pdtc(ig,l)=0. ! final tendancy temperature set to 0
251
252	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
253	condsub(ig)=.true. !condensation at level l
254	IF (zfallice(ig,l+1).gt.0) then
255	zfallheat=zfallice(ig,l+1)*&
256	(pphi(ig,l+1)-pphi(ig,l) +&
257	cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2
258	ELSE
259	zfallheat=0.
260	ENDIF
261	zdtlatent(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
262	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))&
263	ccondzdtlatent(ig,l)- zfallheat
264	! Case when the ice from above sublimes entirely
265	! """""""""""""""""""""""""""""""""""""""""""""""
266	IF ((zfallice(ig,l+1).lt.-zcondicea(ig,l)) &
267	.AND. (zfallice(ig,l+1).gt.0)) THEN
268
269	zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/&
270	(ccond*(pplev(ig,l)-pplev(ig,l+1)))
271	zcondicea(ig,l)= -zfallice(ig,l+1)
272	END IF
273
274	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
275
276	END IF
277
278	ENDDO
279	ENDDO
280	endif ! not fast
281
282	!-----------------------------------------------------------------------
283	! Condensation/sublimation on the ground
284	! (calculation of zcondices and pdtsrfc)
285
286	! Adding subtimesteps : in fast version, pressures too low lead to
287	! instabilities.
288	IF (fast) THEN
289	IF (pplev(1,1).gt.0.3) THEN
290	nsubtimestep= 1
291	ELSE
292	nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1)
293	ENDIF
294	ELSE
295	nsubtimestep= 1 ! if more then kp and p00 have to be calculated
296	! for nofast mode
297	ENDIF ! fast
298	subtimestep=ptimestep/float(nsubtimestep)
299
300	do itsub=1,nsubtimestep
301	! first loop : getting zplev, ztsurf
302	IF (itsub.eq.1) then
303	DO ig=1,klon
304	zplev(ig)=pplev(ig,1)
305	ztsrfhist(ig)=ptsrf(ig) + pdtsrf(ig)*ptimestep
306	ztsrf(ig)=ptsrf(ig) + pdtsrf(ig)*subtimestep !!
307	zpicen2(ig)=picen2(ig)
308	ENDDO
309	ELSE
310	! additional loop :
311	! 1) pressure updated
312	! 2) direct redistribution of pressure over the globe
313	! 3) modification pressure for unstable cases
314	! 4) pressure update to conserve mass
315	! 5) temperature updated with radiative tendancies
316	DO ig=1,klon
317	zplevhist(ig)=zplev(ig)
318	zplevnew(ig)=zplev(ig)+pdpsrf(ig)*subtimestep ! 1)
319	!IF (zplevnew(ig).le.0.0001) then
320	! zplevnew(ig)=0.0001*kp(ig)/p00
321	!ENDIF
322	ENDDO
323	! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code)
324	call planetwide_sumval(zplevnew(:)*cell_area(:)/totarea_planet,globzplevnew)
325	DO ig=1,klon
326	zplev(ig)=kp(ig)*globzplevnew/p00 ! 2)
327	ENDDO
328	DO ig=1,klon ! 3) unstable case condensation and sublimation: zplev=zplevhist
329	IF (((pdpsrf(ig).lt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).le.ztsrf(ig))).or. &
330	((pdpsrf(ig).gt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).ge.ztsrf(ig)))) then
331	zplev(ig)=zplevhist(ig)
332	ENDIF
333	zplevhist(ig)=zplev(ig)
334	ENDDO
335	call planetwide_sumval(zplevhist(:)*cell_area(:)/totarea_planet,zplevhistglob)
336	zplev=zplev*globzplevnew/zplevhistglob ! 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	if (fast) then
360	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
361	/(lw_n2*subtimestep)
362	else
363	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
364	/((lw_n2+cpp(zt(ig,1)-ztcondsol(ig)))subtimestep)
365	endif
366	else ! sublimation
367	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
368	/(lw_n2*subtimestep)
369	end if
370
371	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep
372
373	! partial sublimation of N2 ice
374	! If the entire N_2 ice layer sublimes
375	! (including what has just condensed in the atmosphere)
376	IF((zpicen2(ig)/subtimestep).LE. &
377	-zcondices(ig))THEN
378	zcondices(ig) = -zpicen2(ig)/subtimestep
379	pdtsrfc(ig)=(lw_n2/pcapcal(ig))* &
380	(zcondices(ig))
381	END IF
382
383	! Changing N2 ice amount and pressure
384	pdicen2(ig) = zcondices(ig)
385	pdpsrf(ig) = -pdicen2(ig)*g
386	! pdpsrf(ig) = 0. ! OPTION to check impact N2 sub/cond
387	IF (fast.and.(zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001)) then
388	pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
389	pdicen2(ig)=-pdpsrf(ig)/g
390	ENDIF
391
392	ELSE ! no condsub
393
394	!IF (zpicen2(ig) .GT. 0.) THEN
395	! print,'TB24 no condsub',latitude(ig)180./pi,longitude(ig)*180./pi
396	!ENDIF
397	pdpsrf(ig)=0.
398	pdicen2(ig)=0.
399	pdtsrfc(ig)=0.
400	ENDIF
401	ENDDO ! ig
402	enddo ! itsub,nsubtimestep
403
404	! Updating pressure, temperature and ice reservoir
405	DO ig=1,klon
406	pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep
407	! Two options here : 1 ok, 2 is wrong
408	pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep
409	!pdicen2(ig)=-pdpsrf(ig)/g
410
411	pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep
412
413	!IF (picen2(ig) .GT. 0. .and. abs(pdicen2(ig)).lt.1.e-12) THEN
414	! print,'TB24 low pdq',latitude(ig)180./pi,longitude(ig)*180./pi
415	!ENDIF
416
417
418
419	! security
420	if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then
421	write(,) 'WARNING in condense_n2:'
422	write(,) picen2(ig),pdicen2(ig)*ptimestep
423	pdicen2(ig)= -picen2(ig)/ptimestep
424	pdpsrf(ig)=-pdicen2(ig)*g
425	endif
426
427	if(.not.picen2(ig).ge.0.) THEN
428	! if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
429	print, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)ptimestep
430	! pdicen2(ig)= -picen2(ig)/ptimestep
431	! else
432	picen2(ig)=0.0
433	! endif
434	endif
435	ENDDO ! ig
436
437	! ***************************************************************
438	! Correction to account for redistribution between sigma or hybrid
439	! layers when changing surface pressure (and warming/cooling
440	! of the n2 currently changing phase).
441	! *************************************************************
442	if (.not.fast) then
443	DO ig=1,klon
444	if (condsub(ig)) then
445
446	! Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
447	! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
448	zmflux(1) = -zcondices(ig)
449	DO l=1,klev
450	zmflux(l+1) = zmflux(l) -zcondicea(ig,l) &
451	+ (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
452	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
453	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
454	END DO
455
456	! Mass of each layer
457	! ------------------
458	DO l=1,klev
459	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
460	END DO
461
462	! Corresponding fluxes for T,U,V,Q
463	! """"""""""""""""""""""""""""""""
464	! averaging operator for TRANSPORT
465	! """"""""""""""""""""""""""""""""
466
467	! Subtimestep loop to perform the redistribution gently and simultaneously with
468	! the other tendencies
469	! Estimation of subtimestep (using only the first layer, the most critical)
470	dtmax=ptimestep
471	if (zmflux(1).gt.1.e-20) then
472	dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
473	endif
474	nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
475	subtimestep=ptimestep/float(nsubtimestep)
476
477	! New flux for each subtimestep
478	do l=1,klev+1
479	w(l)=-zmflux(l)*subtimestep
480	enddo
481	! initializing variables that will vary during subtimestep:
482	do l=1,klev
483	ztc(l) =pt(ig,l)
484	zu(l) =pu(ig,l)
485	zv(l) =pv(ig,l)
486	do iq=1,nq
487	zq(l,iq) = pq(ig,l,iq)
488	enddo
489	end do
490
491	! loop over nsubtimestep
492	! """"""""""""""""""""""
493	do itsub=1,nsubtimestep
494	! Progressively adding tendancies from other processes.
495	do l=1,klev
496	ztc(l) =ztc(l) +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
497	zu(l) =zu(l) +pdu( ig,l) * subtimestep
498	zv(l) =zv(l) +pdv( ig,l) * subtimestep
499	do iq=1,nq
500	zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
501	enddo
502	end do
503
504	! Change of mass in each layer
505	do l=1,klev
506	masse(l)=masse(l)+pdpsrf(ig)subtimestep(pplev(ig,l) - pplev(ig,l+1))&
507	/(g*pplev(ig,1))
508	end do
509
510	ztm(:)=0.
511	zum(:)=0.
512	zvm(:)=0.
513	zqm1(:)=0.
514	zqm(:,:)=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	print*,"m", m
824	print*,"Mtot", Mtot
825	print*,"Mqtot", Mqtot
826	call abort_physic("condense_n2","top layer is disappearing !",1)
827	end if
828	else ! if(w(l+1).lt.0)
829	m = l-1
830	Mtot = masse(m+1)
831	MQtot = masse(m+1)*q(m+1)
832	if (m.gt.0) then ! because some compilers will have problems
833	! evaluating masse(0)
834	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
835	m=m-1
836	Mtot = Mtot + masse(m+1)
837	MQtot = MQtot + masse(m+1)*q(m+1)
838	if (m.eq.0) exit
839	end do
840	endif
841	if (m.gt.0) then
842	sigw=(w(l+1)+Mtot)/masse(m)
843	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot) &
844	(q(m)-0.5(1.+sigw)dzq(m)) )
845	else
846	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
847	end if
848	end if
849	enddo
850
851	return
852	end subroutine vl1d

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

Download in other formats:

Original Format