Context Navigation

vlz_fi.F @ 1655

Last change on this file since 1655 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

File size: 5.6 KB

Line
1	SUBROUTINE vlz_fi(ngrid,nlay,q,pente_max,masse,w,wq)
2	c
3	c Auteurs: P.Le Van, F.Hourdin, F.Forget
4	c
5	c ********************************************************************
6	c Shema d'advection " pseudo amont " dans la verticale
7	c pour appel dans la physique (sedimentation)
8	c ********************************************************************
9	c q rapport de melange (kg/kg)...
10	c masse : masse de la couche Dp/g
11	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
12	c pente_max = 2 conseillee
13	c
14	c
15	c --------------------------------------------------------------------
16	IMPLICIT NONE
17	c
18
19	c
20	c
21	c Arguments:
22	c ----------
23	integer,intent(in) :: ngrid ! number of atmospheric columns
24	integer,intent(in) :: nlay ! number of atmospheric layers
25	real masse(ngrid,nlay),pente_max
26	REAL q(ngrid,nlay)
27	REAL w(ngrid,nlay)
28	REAL wq(ngrid,nlay+1)
29	c
30	c Local
31	c ---------
32	c
33	INTEGER i,ij,l,j,ii
34	c
35
36	real dzq(ngrid,nlay),dzqw(ngrid,nlay),adzqw(ngrid,nlay),dzqmax
37	real newmasse
38	real sigw, Mtot, MQtot
39	integer m
40
41	REAL SSUM,CVMGP,CVMGT
42	integer ismax,ismin
43
44
45	c On oriente tout dans le sens de la pression c'est a dire dans le
46	c sens de W
47
48	do l=2,nlay
49	do ij=1,ngrid
50	dzqw(ij,l)=q(ij,l-1)-q(ij,l)
51	adzqw(ij,l)=abs(dzqw(ij,l))
52	enddo
53	enddo
54
55	do l=2,nlay-1
56	do ij=1,ngrid
57	#ifdef CRAY
58	dzq(ij,l)=0.5*
59	, cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
60	#else
61	if(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) then
62	dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
63	else
64	dzq(ij,l)=0.
65	endif
66	#endif
67	dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
68	dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
69	enddo
70	enddo
71
72	do ij=1,ngrid
73	dzq(ij,1)=0.
74	dzq(ij,nlay)=0.
75	enddo
76	c ---------------------------------------------------------------
77	c .... calcul des termes d'advection verticale .......
78	c ---------------------------------------------------------------
79
80	c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
81	c
82	c No flux at the model top:
83	do ij=1,ngrid
84	wq(ij,nlay+1)=0.
85	enddo
86
87	c 1) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION)
88	c ===============================
89
90	do l = 1,nlay ! loop different than when w<0
91	do ij = 1,ngrid
92
93	if(w(ij,l).gt.0.)then
94
95	c Regular scheme (transfered mass < 1 layer)
96	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
97	if(w(ij,l).le.masse(ij,l))then
98	sigw=w(ij,l)/masse(ij,l)
99	wq(ij,l)=w(ij,l)(q(ij,l)+0.5(1.-sigw)*dzq(ij,l))
100
101
102	c Extended scheme (transfered mass > 1 layer)
103	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
104	else
105	m=l
106	Mtot = masse(ij,m)
107	MQtot = masse(ij,m)*q(ij,m)
108	if(m.ge.nlay)goto 88
109	do while(w(ij,l).gt.(Mtot+masse(ij,m+1)))
110	m=m+1
111	Mtot = Mtot + masse(ij,m)
112	MQtot = MQtot + masse(ij,m)*q(ij,m)
113	if(m.ge.nlay)goto 88
114	end do
115	88 continue
116	if (m.lt.nlay) then
117	sigw=(w(ij,l)-Mtot)/masse(ij,m+1)
118	wq(ij,l)=(MQtot + (w(ij,l)-Mtot)*
119	& (q(ij,m+1)+0.5(1.-sigw)dzq(ij,m+1)) )
120	else
121	w(ij,l) = Mtot
122	wq(ij,l) = Mqtot
123	end if
124	end if
125	end if
126	enddo
127	enddo
128
129	c 2) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION)
130	c ===============================
131	goto 99 ! SKIPPING THIS PART FOR SEDIMENTATION
132
133	c Surface flux up:
134	do ij = 1,ngrid
135	if(w(ij,1).lt.0.) wq(ij,1)=0. ! warning : not always valid
136	end do
137
138	do l = 1,nlay-1 ! loop different than when w>0
139	do ij = 1,ngrid
140	if(w(ij,l+1).le.0)then
141
142	c Regular scheme (transfered mass < 1 layer)
143	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
144	if(-w(ij,l+1).le.masse(ij,l))then
145	sigw=w(ij,l+1)/masse(ij,l)
146	wq(ij,l+1)=w(ij,l+1)(q(ij,l)-0.5(1.+sigw)*dzq(ij,l))
147	c Extended scheme (transfered mass > 1 layer)
148	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
149	else
150	m = l-1
151	Mtot = masse(ij,m+1)
152	MQtot = masse(ij,m+1)*q(ij,m+1)
153	if (m.le.0)goto 77
154	do while(-w(ij,l+1).gt.(Mtot+masse(ij,m)))
155	m=m-1
156	Mtot = Mtot + masse(ij,m+1)
157	MQtot = MQtot + masse(ij,m+1)*q(ij,m+1)
158	if (m.le.0)goto 77
159	end do
160	77 continue
161
162	if (m.gt.0) then
163	sigw=(w(ij,l+1)+Mtot)/masse(ij,m)
164	wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)*
165	& (q(ij,m)-0.5(1.+sigw)dzq(ij,m)) )
166	else
167	c wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)*qm(ij,1))
168	write(,) 'a rather weird situation in vlz_fi !'
169	stop
170	end if
171	endif
172	endif
173	enddo
174	enddo
175	99 continue
176
177	do l=1,nlay
178	do ij=1,ngrid
179
180	cccccccc lines below not used for sedimentation (No real flux)
181	ccccc newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
182	ccccc q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
183	ccccc& /newmasse
184	ccccc masse(ij,l)=newmasse
185
186	q(ij,l)=q(ij,l) + (wq(ij,l+1)-wq(ij,l))/masse(ij,l)
187
188	enddo
189	enddo
190
191
192
193	return
194	end

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

Download in other formats:

Original Format