Context Navigation

newcondens.F @ 1635

Last change on this file since 1635 was 1621, checked in by emillour, 8 years ago

Further work on full dynamics/physics separation.

LMDZ.COMMON:

added phy_common/vertical_layers_mod.F90 to store information on vertical grid. This is where routines in the physics should get the information.
The contents of vertical_layers_mod intialized via dynphy_lonlat/inigeomphy_mod.F90.

LMDZ.MARS:

physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
created an "ini_tracer_mod" routine in module "tracer_mod" for a cleaner initialization of the later.
removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.GENERIC:

physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
added nqtot to tracer_h.F90.
removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.VENUS:

physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics. Initialized via iniphysiq. IMPORTANT: there are some hard-coded constants! These should match what is in cpdet_mod.F90 in the dynamics.
got rid of references to moyzon_mod module within the physics. The required variables (tmoy, plevmoy) are passed to the physics as arguments to physiq.

LMDZ.TITAN:

added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics.
Extra work required to completely decouple physics and dynamics: moyzon_mod should be cleaned up and information passed from dynamics to physics as as arguments. Likewise moyzon_ch and moyzon_mu should not be queried from logic_mod (which is in the dynamics).

File size: 33.1 KB

Line
1	SUBROUTINE newcondens(ngrid,nlayer,nq,ptimestep,
2	$ pcapcal,pplay,pplev,ptsrf,pt,
3	$ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
4	$ piceco2,psolaralb,pemisurf,
5	$ pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
6	$ fluxsurf_sw,zls)
7
8	use tracer_mod, only: noms
9	use surfdat_h, only: emissiv, phisfi
10	use geometry_mod, only: latitude ! grid point latitudes (rad)
11	use planete_h
12	USE comcstfi_h
13	#ifndef MESOSCALE
14	USE vertical_layers_mod, ONLY: bp
15	#endif
16	IMPLICIT NONE
17	c=======================================================================
18	c subject:
19	c --------
20	c Condensation/sublimation of CO2 ice on the ground and in the
21	c atmosphere
22	c (Scheme described in Forget et al., Icarus, 1998)
23	c
24	c author: Francois Forget 1994-1996
25	c ------
26	c
27	c input:
28	c -----
29	c ngrid nombre de points de verticales
30	c (toutes les boucles de la physique sont au
31	c moins vectorisees sur ngrid)
32	c nlayer nombre de couches
33	c pplay(ngrid,nlayer) Pressure levels
34	c pplev(ngrid,nlayer+1) Pressure levels
35	c pt(ngrid,nlayer) temperature (en K)
36	c ptsrf(ngrid) temperature de surface
37	c
38	c \
39	c pdt(ngrid,nlayer)\ derivee temporelle physique avant condensation
40	c / ou sublimation pour pt,ptsrf
41	c pdtsrf(ngrid) /
42	c
43	c output:
44	c -------
45	c
46	c pdpsrf(ngrid) \ derivee temporelle physique (contribution de
47	c pdtc(ngrid,nlayer) / la condensation ou sublimation) pour Ps,pt,ptsrf
48	c pdtsrfc(ngrid) /
49	c
50	c Entree/sortie
51	c -------------
52	c
53	c piceco2(ngrid) : quantite de glace co2 au sol (kg/m2)
54	c psolaralb(ngrid,2) : albedo au sol
55	c pemisurf(ngrid) : emissivite du sol
56
57	c
58	c=======================================================================
59	c
60	c 0. Declarations :
61	c ------------------
62	c
63	include "callkeys.h"
64
65	c-----------------------------------------------------------------------
66	c Arguments :
67	c ---------
68	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
69	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
70	INTEGER,INTENT(IN) :: nq ! number of tracers
71
72	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
73	REAL,INTENT(IN) :: pcapcal(ngrid)
74	REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
75	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
76	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
77	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
78	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
79	REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
80	! previous physical processes (K/s)
81	REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
82	! from previous physical processes
83	REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
84	! from previous physical processes
85	REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from
86	! previous physical processes (K/s)
87	REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
88	REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
89	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
90	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
91	! previous physical processes
92	REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2)
93	REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface
94	REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface
95
96	! tendencies due to CO2 condensation/sublimation:
97	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
98	REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s)
99	REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
100	REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
101	REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
102	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
103
104	! added to calculate flux dependent albedo:
105	REAL,intent(in) :: fluxsurf_sw(ngrid,2)
106	real,intent(in) :: zls ! solar longitude (rad)
107
108	c
109	c Local variables :
110	c -----------------
111
112	c variables used for albedo parametrization
113	c --------------------------------------------
114	INTEGER i,j
115	c REAL Fluxmean(jjp1)
116	INTEGER l,ig,iq,icap,nmix
117	LOGICAL transparency, fluxdependent
118	c flag transparency if you want to make the co2ice semi-transparent
119	PARAMETER(transparency=.true.)
120	c flag fluxdependent if you want the co2ice albedo to be dependent on
121	c the incident solar flux
122	PARAMETER(fluxdependent=.false.)
123	REAL slopy,alpha,constA,constB,constT,albediceF_new(ngrid)
124	REAL zt(ngrid,nlayer)
125	REAL zcpi
126	REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm)
127	REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface)
128	REAL zdiceco2(ngrid)
129	REAL zcondicea(ngrid,nlayer)
130	REAL zcondices(ngrid)
131	REAL zfallice(ngrid,nlayer+1) , zfallheat
132	REAL zmflux(nlayer+1)
133	REAL zu(nlayer),zv(nlayer)
134	REAL zq(nlayer,nq),zq1(nlayer)
135	REAL ztsrf(ngrid)
136	REAL ztc(nlayer), ztm(nlayer+1)
137	REAL zum(nlayer+1) , zvm(nlayer+1)
138	REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
139	REAL masse(nlayer),w(nlayer+1)
140	REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
141	LOGICAL condsub(ngrid)
142
143	real :: emisref(ngrid)
144
145	c variable speciale diagnostique
146	real tconda1(ngrid,nlayer)
147	real tconda2(ngrid,nlayer)
148	c REAL zdiceco2a(ngrid) ! for diagnostic only
149	real zdtsig (ngrid,nlayer)
150	real zdt (ngrid,nlayer)
151	real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
152	! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
153	! then condensation temperature is computed using partial pressure of CO2
154	logical,parameter :: improved_ztcond=.true.
155	! Bound co2 (tracer) values...
156	logical,parameter :: bound_qco2=.false.
157	real,parameter :: qco2max=1.1
158	real,parameter :: qco2mini=0.1
159	real :: zqco2
160
161	c local saved variables
162	integer,save :: ico2 ! index of CO2 tracer
163	real,save :: qco2min,qco2,mmean
164	real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
165	real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
166	real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
167	REAL,SAVE :: acond,bcond,ccond
168	! REAL,SAVE :: albediceF(ngrid)
169	real,save :: m_co2, m_noco2, A , B
170
171	LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
172
173	integer flag
174
175	c----------------------------------------------------------------------
176
177	c Initialisation
178	c --------------
179	c
180	IF (firstcall) THEN
181
182	bcond=1./tcond1mb
183	ccond=cpp/(g*latcond)
184	acond=r/latcond
185
186	firstcall=.false.
187	write(,) 'Newcondens: improved_ztcond=',improved_ztcond
188	write(,) 'Newcondens: bound_qco2=',bound_qco2
189	PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
190	PRINT*,'acond,bcond,ccond',acond,bcond,ccond
191
192	ico2=0
193
194	if (tracer) then
195	c Prepare Special treatment if one of the tracer is CO2 gas
196	do iq=1,nq
197	if (noms(iq).eq."co2") then
198	ico2=iq
199	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
200	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
201	c Compute A and B coefficient use to compute
202	c mean molecular mass Mair defined by
203	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
204	c 1/Mair = A*q(ico2) + B
205	A =(1/m_co2 - 1/m_noco2)
206	B=1/m_noco2
207	endif
208	enddo
209	c minimum CO2 mix. ratio below which mixing occur with layer above:
210	qco2min =0.75
211	end if
212	ENDIF ! of IF (firstcall)
213	zcpi=1./cpp
214
215	c
216	c======================================================================
217	c Calcul of CO2 condensation sublimation
218	c ============================================================
219	c
220	c Used variable :
221	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
222	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
223	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
224	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
225	c
226	c pdtc(ngrid,nlayer) : dT/dt due to cond/sub
227	c
228	c
229	c Tendencies set to 0 (except pdtc)
230	c -------------------------------------
231	DO l=1,nlayer
232	DO ig=1,ngrid
233	zcondicea(ig,l) = 0.
234	zfallice(ig,l) = 0.
235	pduc(ig,l) = 0
236	pdvc(ig,l) = 0
237	END DO
238	END DO
239
240	DO iq=1,nq
241	DO l=1,nlayer
242	DO ig=1,ngrid
243	pdqc(ig,l,iq) = 0
244	END DO
245	END DO
246	END DO
247
248	DO ig=1,ngrid
249	zfallice(ig,nlayer+1) = 0.
250	zcondices(ig) = 0.
251	pdtsrfc(ig) = 0.
252	pdpsrf(ig) = 0.
253	condsub(ig) = .false.
254	zdiceco2(ig) = 0.
255	ENDDO
256	zfallheat=0
257
258	c *************************
259	c ATMOSPHERIC CONDENSATION
260	c *************************
261
262	c Compute CO2 Volume mixing ratio
263	c -------------------------------
264	if (improved_ztcond.and.(ico2.ne.0)) then
265	DO l=1,nlayer
266	DO ig=1,ngrid
267	qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
268	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
269	mmean=1/(A*qco2 +B)
270	vmr_co2(ig,l) = qco2*mmean/m_co2
271	ENDDO
272	ENDDO
273	else
274	DO l=1,nlayer
275	DO ig=1,ngrid
276	vmr_co2(ig,l)=0.95
277	ENDDO
278	ENDDO
279	end if
280
281	c forecast of atmospheric temperature zt and frost temperature ztcond
282	c --------------------------------------------------------------------
283
284	DO l=1,nlayer
285	DO ig=1,ngrid
286	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
287	! ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
288	if (pplay(ig,l).ge.1e-4) then
289	ztcond(ig,l)=
290	& 1./(bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
291	else
292	ztcond(ig,l)=0.0 !mars Monica
293	endif
294	ENDDO
295	ENDDO
296
297	ztcond(:,nlayer+1)=ztcond(:,nlayer)
298
299	c Condensation/sublimation in the atmosphere
300	c ------------------------------------------
301	c (calcul of zcondicea , zfallice and pdtc)
302	c
303	DO l=nlayer , 1, -1
304	DO ig=1,ngrid
305	pdtc(ig,l)=0.
306	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
307	condsub(ig)=.true.
308	IF (zfallice(ig,l+1).gt.0) then
309	zfallheat=zfallice(ig,l+1)*
310	& (pphi(ig,l+1)-pphi(ig,l) +
311	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
312	ELSE
313	zfallheat=0.
314	ENDIF
315	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
316	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
317	& ccondpdtc(ig,l)- zfallheat
318	c Case when the ice from above sublimes entirely
319	c """""""""""""""""""""""""""""""""""""""""""""""
320	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
321	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
322	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
323	zcondicea(ig,l)= -zfallice(ig,l+1)
324	END IF
325
326	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
327	END IF
328	ENDDO
329	ENDDO
330
331	c *************************
332	c SURFACE CONDENSATION
333	c *************************
334
335	c forecast of ground temperature ztsrf and frost temperature ztcondsol
336	c --------------------------------------------------------------------
337	DO ig=1,ngrid
338	ztcondsol(ig)=
339	& 1./(bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
340	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
341	ENDDO
342
343	c
344	c Condensation/sublimation on the ground
345	c --------------------------------------
346	c (compute zcondices and pdtsrfc)
347	c
348	DO ig=1,ngrid
349	IF(latitude(ig).lt.0) THEN
350	! Southern hemisphere
351	icap=2
352	ELSE
353	! Northern hemisphere
354	icap=1
355	ENDIF
356
357	c
358	c Loop on where we have condensation/ sublimation
359	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
360	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
361	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
362	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
363	condsub(ig) = .true.
364
365	IF (zfallice(ig,1).gt.0) then
366	zfallheat=zfallice(ig,1)*
367	& (pphi(ig,1)- phisfi(ig) +
368	& cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
369	ELSE
370	zfallheat=0.
371	ENDIF
372
373	c condensation or partial sublimation of CO2 ice
374	c """""""""""""""""""""""""""""""""""""""""""""""
375	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
376	& /(latcond*ptimestep) - zfallheat
377	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
378
379	c If the entire CO_2 ice layer sublimes
380	c """"""""""""""""""""""""""""""""""""""""""""""""""""
381	c (including what has just condensed in the atmosphere)
382
383	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
384	& -zcondices(ig))THEN
385	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
386	pdtsrfc(ig)=(latcond/pcapcal(ig))*
387	& (zcondices(ig)+zfallheat)
388	END IF
389
390	c Changing CO2 ice amount and pressure :
391	c """"""""""""""""""""""""""""""""""""
392
393	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
394	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
395	pdpsrf(ig) = -zdiceco2(ig)*g
396
397	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
398	PRINT*,'STOP in condens'
399	PRINT*,'condensing more than total mass'
400	PRINT*,'Grid point ',ig
401	PRINT*,'Ps = ',pplev(ig,1)
402	PRINT*,'d Ps = ',pdpsrf(ig)
403	STOP
404	ENDIF
405	END IF ! if there is condensation/sublimmation
406	ENDDO ! of DO ig=1,ngrid
407
408	c ********************************************************************
409	c Surface albedo and emissivity of the surface below the snow (emisref)
410	c ********************************************************************
411	c Prepare the case where albedo varies with insolation:
412	c ----------------------------------------------------
413	! if (fluxdependent) then
414	!
415	c Calcul du flux moyen (zonal mean)
416	! do j=1,jjp1
417	! Fluxmean(j)=0
418	! do i=1,iim
419	! ig=1+(j-2)*iim +i
420	! if(j.eq.1) ig=1
421	! if(j.eq.jjp1) ig=ngrid
422	! Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
423	! $ +fluxsurf_sw(ig,2)
424	! enddo
425	! Fluxmean(j)=Fluxmean(j)/float(iim)
426	! enddo
427	!
428	c const A and B used to calculate the albedo which depends on solar flux
429	c albedice=constA+constB*Flux
430	c constT = time step to calculate the solar flux when flux decreases
431	! constA=0.26
432	c constA=0.33
433	c constA=0.186
434	! constB=0.00187
435	! constT=10
436	! endif ! of if (fluxdependent)
437
438	! Check that amont of CO2 ice is not problematic
439	DO ig=1,ngrid
440	if(.not.piceco2(ig).ge.0.) THEN
441	if(piceco2(ig).le.-5.e-8) print*,
442	$ 'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
443	piceco2(ig)=0.
444	endif
445	ENDDO
446
447	! Set albedo and emissivity of the surface
448	! ----------------------------------------
449	CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
450
451	c Calcul de l'albedo
452	c ------------------
453	! do ig =1,ngrid
454	! IF(ig.GT.ngrid/2+1) THEN
455	! icap=2
456	! ELSE
457	! icap=1
458	! ENDIF
459	! IF(firstcall2) THEN
460	! albediceF(ig)=albedice(icap)
461	! ENDIF
462	c if there is still co2ice ccccccccccccccccccccccc
463	! if (piceco2(ig).gt.0) then
464	! emisref(ig) = emisice(icap)
465
466	c if flux dependent albedo is used
467	c --------------------------------
468	! if (fluxdependent) then
469	! j=INT((ig-2)/iim)+2
470	! if(ig.eq.1) j=1
471	! if(ig.eq.ngrid) j=jjp1
472	c albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
473	c $ constA+constB*250)
474	! albediceF_new(ig)=constA+constB*Fluxmean(j)
475	! if (albediceF(ig).gt.albediceF_new(ig)) then
476	! albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
477	! $ constT)*(albediceF_new(ig)-albediceF(ig))
478	! else
479	! albediceF(ig)=albediceF_new(ig)
480	! endif
481	c if part of the ice is transparent
482	c slopy = pente de la droite: alpha = y*co2ice/1620
483	c pour une valeur superieur a une epaisseur de glace donnee
484	c ici, epaisseur limite = 10cm
485	! if (transparency) then
486	! slopy=1/(1620*0.10)
487	! alpha=MIN(slopy*piceco2(ig),1.)
488	! psolaralb(ig,1) = alpha*albediceF(ig)
489	! $ +(1-alpha)*albedodat(ig)
490	! psolaralb(ig,2) = psolaralb(ig,1)
491	! else
492	! psolaralb(ig,1) = albediceF(ig)
493	! psolaralb(ig,2) = psolaralb(ig,1)
494	! endif
495	! else
496	c transparency set to true and fluxdependent set to false
497	! if (transparency) then
498	! slopy=1/(1620*0.10)
499	! alpha=MIN(slopy*piceco2(ig),1.)
500	! psolaralb(ig,1) = alpha*albedice(icap)
501	! $ +(1-alpha)*albedodat(ig)
502	! psolaralb(ig,2) = psolaralb(ig,1)
503	! else
504	c simplest case: transparency and flux dependent set to false
505	! psolaralb(ig,1) = albedice(icap)
506	! psolaralb(ig,2) = albedice(icap)
507	! endif
508	! endif
509	c no more co2ice, albedo = ground albedo
510	! else
511	! psolaralb(ig,1) = albedodat(ig)
512	! psolaralb(ig,2) = albedodat(ig)
513	! emisref(ig) = emissiv
514	! pemisurf(ig) = emissiv
515	! endif
516	! end do ! end of the ig loop
517
518	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
519	DO ig=1,ngrid
520	if (piceco2(ig).eq.0) then
521	pemisurf(ig)=emissiv
522	endif
523	ENDDO
524
525	! firstcall2=.false.
526	c ***************************************************************
527	c Correction to account for redistribution between sigma or hybrid
528	c layers when changing surface pressure (and warming/cooling
529	c of the CO2 currently changing phase).
530	c *************************************************************
531
532	DO ig=1,ngrid
533	if (condsub(ig)) then
534	do l=1,nlayer
535	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
536	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
537	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
538	do iq=1,nq
539	zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
540	enddo
541	end do
542
543	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
544	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
545
546	zmflux(1) = -zcondices(ig)
547	DO l=1,nlayer
548	zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
549	#ifndef MESOSCALE
550	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
551	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
552	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
553	#else
554	if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
555	#endif
556	END DO
557
558	c Mass of each layer
559	c ------------------
560	DO l=1,nlayer
561	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
562	END DO
563
564
565	c Corresponding fluxes for T,U,V,Q
566	c """"""""""""""""""""""""""""""""
567
568	c averaging operator for TRANSPORT
569	c """"""""""""""""""""""""""""""""
570	c Value transfert at the surface interface when condensation
571	c sublimation:
572	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
573	zum(1) = 0
574	zvm(1) = 0
575	do iq=1,nq
576	zqm(1,iq)=0. ! most tracer do not condense !
577	enddo
578	c Special case if one of the tracer is CO2 gas
579	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
580
581	c Van Leer scheme:
582	DO l=1,nlayer+1
583	w(l)=-zmflux(l)*ptimestep
584	END DO
585	call vl1d(nlayer,ztc,2.,masse,w,ztm)
586	call vl1d(nlayer,zu ,2.,masse,w,zum)
587	call vl1d(nlayer,zv ,2.,masse,w,zvm)
588	do iq=1,nq
589	do l=1,nlayer
590	zq1(l)=zq(l,iq)
591	enddo
592	zqm1(1)=zqm(1,iq)
593	call vl1d(nlayer,zq1,2.,masse,w,zqm1)
594	do l=2,nlayer
595	zq( l,iq)=zq1(l)
596	zqm(l,iq)=zqm1(l)
597	enddo
598	enddo
599
600	c Surface condensation affects low winds
601	if (zmflux(1).lt.0) then
602	zum(1)= zu(1) * (w(1)/masse(1))
603	zvm(1)= zv(1) * (w(1)/masse(1))
604	if (w(1).gt.masse(1)) then ! ensure numerical stability
605	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
606	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
607	end if
608	end if
609
610	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
611	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
612	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
613	do iq=1,nq
614	zqm(nlayer+1,iq)= zq(nlayer,iq)
615	enddo
616
617	#ifdef MESOSCALE
618	!!!! AS: This part must be commented in the mesoscale model
619	!!!! AS: ... to avoid instabilities.
620	!!!! AS: you have to compile with -DMESOSCALE to do so
621	#else
622	c Tendencies on T, U, V, Q
623	c """"""""""""""""""""""""
624	DO l=1,nlayer
625
626	c Tendencies on T
627	zdtsig(ig,l) = (1/masse(l)) *
628	& ( zmflux(l)*(ztm(l) - ztc(l))
629	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
630	& + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
631	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
632
633	c Tendencies on U
634	pduc(ig,l) = (1/masse(l)) *
635	& ( zmflux(l)*(zum(l) - zu(l))
636	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
637
638
639	c Tendencies on V
640	pdvc(ig,l) = (1/masse(l)) *
641	& ( zmflux(l)*(zvm(l) - zv(l))
642	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
643
644	END DO
645
646	#endif
647
648	c Tendencies on Q
649	do iq=1,nq
650	! if (noms(iq).eq.'co2') then
651	if (iq.eq.ico2) then
652	c SPECIAL Case when the tracer IS CO2 :
653	DO l=1,nlayer
654	pdqc(ig,l,iq)= (1/masse(l)) *
655	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
656	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
657	& + zcondicea(ig,l)*(zq(l,iq)-1.) )
658	END DO
659	else
660	DO l=1,nlayer
661	pdqc(ig,l,iq)= (1/masse(l)) *
662	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
663	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
664	& + zcondicea(ig,l)*zq(l,iq) )
665	END DO
666	end if
667	enddo
668
669	c --------------------------------------------------------
670	c Roughly Simulate Molecular mixing when CO2 is too depleted by
671	c Surface condensation (mixing starts if qco2 < qco2min )
672	c FF 06/2004
673	c WARNING : this is now done in convadj, better (FF 02/2005)
674	c --------------------------------------------------------
675	flag=0 ! now done in convadj : must be =0
676	if (flag.eq.1) then
677	if(ico2.gt.0) then ! relevant only if one tracer is CO2
678	if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
679	& .lt.qco2min) then
680	do iq=1,nq
681	zq(1,iq)=pq(ig,1,iq)
682	& +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
683	Smq(1,iq) = masse(1)*zq(1,iq)
684	end do
685	Sm(1) = masse(1)
686	do l =2,nlayer
687	do iq=1,nq
688	zq(l,iq)=pq(ig,l,iq)
689	& +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
690	smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
691	end do
692	sm(l) = sm(l-1) + masse(l)
693	if(zq(l,ico2).gt.qco2min) then
694	c mixmas: mass of atmosphere that must be mixed to reach qco2min
695	mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
696	& /(zq(l,ico2)-qco2min)
697	if((mixmas.le.sm(l)))then
698	c OK if mixed mass less than mass of layers affected
699	nmix=l ! number of layer affected by mixing
700	goto 99
701	end if
702	end if
703	end do
704	99 continue
705	do iq=1,nq
706	qmix=zq(nmix,iq)
707	& +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
708	do l=1,nmix-1
709	pdqc(ig,l,iq)=
710	& (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
711	end do
712	c layer only partly mixed :
713	pdqc(ig,nmix,iq)=(
714	& qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
715	& -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
716
717	end do
718	end if
719	end if
720
721	endif ! (flag.eq.1)
722	end if ! if (condsub)
723	END DO ! loop on ig
724
725	c ***************************************************************
726	c CO2 snow / clouds scheme
727	c ***************************************************************
728
729	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
730	& zcondicea,zcondices,zfallice,pemisurf)
731
732	c ***************************************************************
733	c Ecriture des diagnostiques
734	c ***************************************************************
735
736	c DO l=1,nlayer
737	c DO ig=1,ngrid
738	c Taux de cond en kg.m-2.pa-1.s-1
739	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
740	c Taux de cond en kg.m-3.s-1
741	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
742	c END DO
743	c END DO
744	c call WRITEDIAGFI(ngrid,'tconda1',
745	c &'Taux de condensation CO2 atmospherique /Pa',
746	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
747	c call WRITEDIAGFI(ngrid,'tconda2',
748	c &'Taux de condensation CO2 atmospherique /m',
749	c & 'kg.m-3.s-1',3,tconda2)
750
751	! output falling co2 ice in 1st layer:
752	! call WRITEDIAGFI(ngrid,'fallice',
753	! &'Precipitation of co2 ice',
754	! & 'kg.m-2.s-1',2,zfallice(1,1))
755
756	!! Specific stuff to bound co2 tracer values ....
757	if (bound_qco2.and.(ico2.ne.0)) then
758	do ig=1,ngrid
759	do l=1,nlayer
760	zqco2=pq(ig,l,ico2)
761	& +(pdq(ig,l,ico2)+pdqc(ig,l,ico2))*ptimestep
762	if (zqco2.gt.qco2max) then
763	! correct pdqc:
764	pdqc(ig,l,ico2)=((qco2max-pq(ig,l,ico2))/ptimestep)
765	& -pdq(ig,l,ico2)
766	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
767	& " so that co2 conc. does not exceed",qco2max
768	write(,) " ig:",ig," l:",l
769	endif ! of if (zqco2.gt.qco2max)
770	if (zqco2.lt.qco2mini) then
771	! correct pdqc:
772	pdqc(ig,l,ico2)=((qco2mini-pq(ig,l,ico2))/ptimestep)
773	& -pdq(ig,l,ico2)
774	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
775	& " so that co2 conc. is not less than",qco2mini
776	write(,) " ig:",ig," l:",l
777	endif ! of if (zqco2.lt.qco2mini)
778	end do
779	enddo
780	endif ! of if (bound_qco2.and.(ico2.ne.0)) then
781
782	#ifndef MESOSCALE
783	! Extra special case for surface temperature tendency pdtsrfc:
784	! we want to fix the south pole temperature to CO2 condensation temperature
785	if (caps.and.(obliquit.lt.27.)) then
786	! check if last grid point is the south pole
787	if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
788	! NB: Updated surface pressure, at grid point 'ngrid', is
789	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
790	! write(,) "newcondens: South pole: latitude(ngrid)=",
791	! & latitude(ngrid)
792	ztcondsol(ngrid)=
793	& 1./(bcond-acondlog(.01vmr_co2(ngrid,1)*
794	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
795	pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
796	endif
797	endif
798	#endif
799
800	return
801	end
802
803
804
805	c *****************************************************************
806	SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm)
807	c
808	c
809	c Operateur de moyenne inter-couche pour calcul de transport type
810	c Van-Leer " pseudo amont " dans la verticale
811	c q,w sont des arguments d'entree pour le s-pg ....
812	c masse : masse de la couche Dp/g
813	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
814	c pente_max = 2 conseillee
815	c
816	c
817	c --------------------------------------------------------------------
818	IMPLICIT NONE
819
820	c
821	c
822	c
823	c Arguments:
824	c ----------
825	integer nlayer
826	real masse(nlayer),pente_max
827	REAL q(nlayer),qm(nlayer+1)
828	REAL w(nlayer+1)
829	c
830	c Local
831	c ---------
832	c
833	INTEGER l
834	c
835	real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
836	real sigw, Mtot, MQtot
837	integer m
838	c integer ismax,ismin
839
840
841	c On oriente tout dans le sens de la pression
842	c W > 0 WHEN DOWN !!!!!!!!!!!!!
843
844	do l=2,nlayer
845	dzqw(l)=q(l-1)-q(l)
846	adzqw(l)=abs(dzqw(l))
847	enddo
848
849	do l=2,nlayer-1
850	if(dzqw(l)*dzqw(l+1).gt.0.) then
851	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
852	else
853	dzq(l)=0.
854	endif
855	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
856	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
857	enddo
858
859	dzq(1)=0.
860	dzq(nlayer)=0.
861
862	do l = 1,nlayer-1
863
864	c Regular scheme (transfered mass < layer mass)
865	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
866	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
867	sigw=w(l+1)/masse(l+1)
868	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
869	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
870	sigw=w(l+1)/masse(l)
871	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
872
873	c Extended scheme (transfered mass > layer mass)
874	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
875	else if(w(l+1).gt.0.) then
876	m=l+1
877	Mtot = masse(m)
878	MQtot = masse(m)*q(m)
879	do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
880	m=m+1
881	Mtot = Mtot + masse(m)
882	MQtot = MQtot + masse(m)*q(m)
883	end do
884	if (m.lt.nlayer) then
885	sigw=(w(l+1)-Mtot)/masse(m+1)
886	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
887	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
888	else
889	w(l+1) = Mtot
890	qm(l+1) = Mqtot / Mtot
891	write(,) 'top layer is disapearing !'
892	stop
893	end if
894	else ! if(w(l+1).lt.0)
895	m = l-1
896	Mtot = masse(m+1)
897	MQtot = masse(m+1)*q(m+1)
898	if (m.gt.0) then ! because some compilers will have problems
899	! evaluating masse(0)
900	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
901	m=m-1
902	Mtot = Mtot + masse(m+1)
903	MQtot = MQtot + masse(m+1)*q(m+1)
904	if (m.eq.0) exit
905	end do
906	endif
907	if (m.gt.0) then
908	sigw=(w(l+1)+Mtot)/masse(m)
909	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
910	& (q(m)-0.5(1.+sigw)dzq(m)) )
911	else
912	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
913	end if
914	end if
915	enddo
916
917	c boundary conditions (not used in newcondens !!)
918	c qm(nlayer+1)=0.
919	c if(w(1).gt.0.) then
920	c qm(1)=q(1)
921	c else
922	c qm(1)=0.
923	c end if
924
925	return
926	end

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

Download in other formats:

Original Format