Context Navigation

condense_n2.F90 @ 3607

Last change on this file since 3607 was 3585, checked in by debatzbr, 2 weeks ago
Connecting microphysics to radiative transfer + miscellaneous cleans
Property svn:executable set to ``*
File size: 27.6 KB

Line
1	subroutine condense_n2(klon,klev,nq,ptimestep, &
2	pcapcal,pplay,pplev,ptsrf,pt, &
3	pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
4	picen2,psolaralb,pemisurf, &
5	pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
6	pdqc,pdicen2)
7
8	use 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 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
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(2:klev+1)=0.
511	zum(2:klev+1)=0.
512	zvm(2:klev+1)=0.
513	zqm1(1:klev+1)=0.
514
515	! Van Leer scheme:
516	call vl1d(klev,ztc,2.,masse,w,ztm)
517	call vl1d(klev,zu ,2.,masse,w,zum)
518	call vl1d(klev,zv ,2.,masse,w,zvm)
519	do iq=1,nq
520	do l=1,klev
521	zq1(l)=zq(l,iq)
522	enddo
523	zqm1(1)=zqm(1,iq)
524	call vl1d(klev,zq1,2.,masse,w,zqm1)
525	do l=2,klev
526	zqm(l,iq)=zqm1(l)
527	enddo
528	enddo
529
530	! Correction to prevent negative value for qn2
531	if (igcm_n2.ne.0) then
532	zqm(1,igcm_n2)=1.
533	do l=1,klev-1
534	if (w(l)zqm(l,igcm_n2).gt.zq(l,igcm_n2)masse(l)) then
535	zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2), &
536	(zqm(l,igcm_n2)w(l) -zq(l,igcm_n2)masse(l))/w(l+1) )
537	else
538	exit
539	endif
540	end do
541	end if
542
543	! Value transfert at the surface interface when condensation sublimation:
544	if (zmflux(1).lt.0) then
545	! Surface condensation
546	zum(1)= zu(1)
547	zvm(1)= zv(1)
548	ztm(1) = ztc(1)
549	else
550	! Surface sublimation:
551	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
552	zum(1) = 0
553	zvm(1) = 0
554	end if
555	do iq=1,nq
556	zqm(1,iq)=0. ! most tracer do not condense !
557	enddo
558	! Special case if the tracer is n2 gas
559	if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
560
561	!!! Source haze: 0.02 pourcent when n2 sublimes
562	IF (source_haze) THEN
563	IF (pdicen2(ig).lt.0) THEN
564	DO iq=1,nq
565	tname=noms(iq)
566	if (tname(1:4).eq."haze") then
567	!zqm(1,iq)=0.02
568	!zqm(1,iq)=-pdicen2(ig)*0.02
569	zqm(1,iq)=-pdicen2(ig)ptimestep0.02
570	!zqm(10,iq)=-pdicen2(ig)ptimestep100.
571	!zqm(1,iq)=-pdicen2(ig)*1000000.
572
573	endif
574	ENDDO
575	ENDIF
576	ENDIF
577	ztm(klev+1)= ztc(klev) ! should not be used, but...
578	zum(klev+1)= zu(klev) ! should not be used, but...
579	zvm(klev+1)= zv(klev) ! should not be used, but...
580	do iq=1,nq
581	zqm(klev+1,iq)= zq(klev,iq)
582	enddo
583
584	! Tendencies on T, U, V, Q
585	! """""""""""""""""""""""
586	DO l=1,klev
587
588	! Tendencies on T
589	zdtsig(ig,l) = (1/masse(l)) * &
590	( zmflux(l)*(ztm(l) - ztc(l)) &
591	- zmflux(l+1)*(ztm(l+1) - ztc(l)) &
592	+ zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
593
594	! Tendencies on U
595	pduc(ig,l) = (1/masse(l)) * &
596	( zmflux(l)*(zum(l) - zu(l))&
597	- zmflux(l+1)*(zum(l+1) - zu(l)) )
598
599	! Tendencies on V
600	pdvc(ig,l) = (1/masse(l)) * &
601	( zmflux(l)*(zvm(l) - zv(l)) &
602	- zmflux(l+1)*(zvm(l+1) - zv(l)) )
603
604	END DO
605
606	! Tendencies on Q
607	do iq=1,nq
608	if (iq.eq.igcm_n2) then
609	! SPECIAL Case when the tracer IS N2 :
610	DO l=1,klev
611	pdqc(ig,l,iq)= (1/masse(l)) * &
612	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
613	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
614	+ zcondicea(ig,l)*(zq(l,iq)-1.) )
615	END DO
616	else
617	DO l=1,klev
618	pdqc(ig,l,iq)= (1/masse(l)) * &
619	( zmflux(l)*(zqm(l,iq) - zq(l,iq)) &
620	- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) &
621	+ zcondicea(ig,l)*zq(l,iq) )
622	END DO
623	end if
624	enddo
625	! Update variables at the end of each subtimestep.
626	do l=1,klev
627	ztc(l) =ztc(l) + zdtsig(ig,l) *subtimestep
628	zu(l) =zu(l) + pduc(ig,l) *subtimestep
629	zv(l) =zv(l) + pdvc(ig,l) *subtimestep
630	do iq=1,nq
631	zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
632	enddo
633	end do
634
635	enddo ! loop on nsubtimestep
636
637	! Recomputing Total tendencies
638	do l=1,klev
639	pdtc(ig,l) = (ztc(l) - zt(ig,l) )/ptimestep
640	pduc(ig,l) = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
641	pdvc(ig,l) = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
642	do iq=1,nq
643	pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep
644
645
646	! Correction temporaire
647	if (iq.eq.igcm_n2) then
648	if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
649	.lt.0.01) then ! if n2 < 1 % !
650	write(,) 'Warning: n2 < 1%'
651	pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
652	end if
653	end if
654
655	enddo
656	end do
657	! *****************************TEMPORAIRE ****************
658	if (klon.eq.1) then
659	write(,) 'nsubtimestep=' ,nsubtimestep
660	write(,) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
661	write(,) 'masse apres' , masse(1)
662	write(,) 'zmfluxDT, l=1' , zmflux(1)ptimestep
663	write(,) 'zmfluxDT, l=2' , zmflux(2)ptimestep
664	write(,) 'pq, l=1,2,3' , pq(1,1,1), pq(1,2,1),pq(1,3,1)
665	write(,) 'zq, l=1,2,3' , zq(1,1), zq(2,1),zq(3,1)
666	write(,) 'dqDt l=1' , pdq(1,1,1)ptimestep
667	write(,) 'dqcondDt l=1' , pdqc(1,1,1)ptimestep
668	end if
669
670	! ***********************************************************
671	end if ! if (condsub)
672	END DO ! loop on ig
673	endif ! not fast
674
675	! ***************************************************************
676	! Ecriture des diagnostiques
677	! ***************************************************************
678
679	! DO l=1,klev
680	! DO ig=1,klon
681	! Taux de cond en kg.m-2.pa-1.s-1
682	! tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
683	! Taux de cond en kg.m-3.s-1
684	! tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
685	! END DO
686	! END DO
687	! call WRITEDIAGFI(klon,'tconda1', &
688	! 'Taux de condensation N2 atmospherique /Pa', &
689	! 'kg.m-2.Pa-1.s-1',3,tconda1)
690	! call WRITEDIAGFI(klon,'tconda2', &
691	! 'Taux de condensation N2 atmospherique /m', &
692	! 'kg.m-3.s-1',3,tconda2)
693
694
695	return
696	end subroutine condense_n2
697
698	!-------------------------------------------------------------------------
699
700	real function tcond_n2(p,vmr)
701	! Calculates the condensation temperature for N2 at pressure P and vmr
702	implicit none
703	real, intent(in):: p,vmr
704
705	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
706	IF (p.ge.0.529995) then
707	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
708	! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e50.98)log(1./(0.1255701.e5)p*vmr))
709	tcond_n2 = (1.)/(0.01583606505-1.211938776e-3log(7.963685594e-5p*vmr))
710	ELSE
711	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
712	! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e51.09)log(1./(0.508059)pvmr))
713	tcond_n2 = (1.)/(1./35.6-1.089633028e-3log(1.968275338p*vmr))
714	ENDIF
715	return
716	end function tcond_n2
717
718	!-------------------------------------------------------------------------
719
720	real function pcond_n2(t,vmr)
721	! Calculates the condensation pressure for N2 at temperature T and vmr
722	implicit none
723	real, intent(in):: t,vmr
724
725	! tcond_n2 = (1.)/(0.026315-0.0011877log(.7143p*vmr))
726	IF (t.ge.35.6) then
727	! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
728	! pcond_n2 = 0.1255701.e5/vmrexp((2.5e50.98)/296.925(1./63.1470-1./t))
729	pcond_n2 = 0.125570e5/vmrexp(825.1241896(1./63.147-1./t))
730	ELSE
731	! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB
732	! pcond_n2 = 0.508059/vmrexp((2.5e51.09)/296.925*(1./35.6-1./t))
733	pcond_n2 = 0.508059/vmrexp(917.7401701(1./35.6-1./t))
734	ENDIF
735	return
736	end function pcond_n2
737
738	!-------------------------------------------------------------------------
739
740	subroutine vl1d(klev,q,pente_max,masse,w,qm)
741	!
742	! Operateur de moyenne inter-couche pour calcul de transport type
743	! Van-Leer " pseudo amont " dans la verticale
744	! q,w sont des arguments d'entree pour le s-pg ....
745	! masse : masse de la couche Dp/g
746	! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
747	! pente_max = 2 conseillee
748	! --------------------------------------------------------------------
749	IMPLICIT NONE
750
751	! Arguments:
752	! ----------
753	integer klev
754	real masse(klev),pente_max
755	REAL q(klev),qm(klev+1)
756	REAL w(klev+1)
757	!
758	! Local
759	! ---------
760	!
761	INTEGER l
762	!
763	real dzq(klev),dzqw(klev),adzqw(klev),dzqmax
764	real sigw, Mtot, MQtot
765	integer m
766
767
768	! On oriente tout dans le sens de la pression
769	! W > 0 WHEN DOWN !!!!!!!!!!!!!
770
771	do l=2,klev
772	dzqw(l)=q(l-1)-q(l)
773	adzqw(l)=abs(dzqw(l))
774	enddo
775
776	do l=2,klev-1
777	if(dzqw(l)*dzqw(l+1).gt.0.) then
778	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
779	else
780	dzq(l)=0.
781	endif
782	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
783	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
784	enddo
785
786	dzq(1)=0.
787	dzq(klev)=0.
788
789	do l = 1,klev-1
790
791	! Regular scheme (transfered mass < layer mass)
792	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
793	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
794	sigw=w(l+1)/masse(l+1)
795	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
796	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
797	sigw=w(l+1)/masse(l)
798	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
799
800	! Extended scheme (transfered mass > layer mass)
801	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
802	else if(w(l+1).gt.0.) then
803	m=l+1
804	Mtot = masse(m)
805	MQtot = masse(m)*q(m)
806	if (m.lt.klev) then ! because some compilers will have problems
807	! evaluating masse(klev+1)
808	do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1))))
809	m=m+1
810	Mtot = Mtot + masse(m)
811	MQtot = MQtot + masse(m)*q(m)
812	if (m.eq.klev) exit
813	end do
814	endif
815	if (m.lt.klev) then
816	sigw=(w(l+1)-Mtot)/masse(m+1)
817	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot) &
818	(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
819	else
820	w(l+1) = Mtot
821	qm(l+1) = Mqtot / Mtot
822	write(,) 'top layer is disapearing !'
823	stop
824	end if
825	else ! if(w(l+1).lt.0)
826	m = l-1
827	Mtot = masse(m+1)
828	MQtot = masse(m+1)*q(m+1)
829	if (m.gt.0) then ! because some compilers will have problems
830	! evaluating masse(0)
831	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
832	m=m-1
833	Mtot = Mtot + masse(m+1)
834	MQtot = MQtot + masse(m+1)*q(m+1)
835	if (m.eq.0) exit
836	end do
837	endif
838	if (m.gt.0) then
839	sigw=(w(l+1)+Mtot)/masse(m)
840	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot) &
841	(q(m)-0.5(1.+sigw)dzq(m)) )
842	else
843	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
844	end if
845	end if
846	enddo
847
848	return
849	end subroutine vl1d

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

Download in other formats:

Original Format