/[lmdze]/trunk/Sources/dyn3d/Vlsplt/vly.f
ViewVC logotype

Contents of /trunk/Sources/dyn3d/Vlsplt/vly.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 32 - (show annotations)
Tue Apr 6 17:52:58 2010 UTC (14 years, 1 month ago) by guez
Original Path: trunk/libf/dyn3d/Vlsplt/vly.f90
File size: 7190 byte(s) 
Split "stringop.f90" into single-procedure files. Gathered files in directory
"IOIPSL/Stringop".

Split "flincom.f90" into "flincom.f90" and "flinget.f90". Removed
unused procedures from module "flincom". Removed unused argument
"filename" of procedure "flinopen_nozoom".

Removed unused files.

Split "grid_change.f90" into "grid_change.f90" and
"gr_phy_write_3d.f90".

Removed unused procedures from modules "calendar", "ioipslmpp",
"grid_atob", "gath_cpl" and "getincom". Removed unused procedures in
files "ppm3d.f" and "thermcell.f".

Split "mathelp.f90" into "mathelp.f90" and "mathop.f90".

Removed unused variable "dpres" of module "comvert".

Use argument "itau" instead of local variables "iadvtr" and "first" to
control algorithm in procedure "fluxstokenc".

Removed unused arguments of procedure "integrd".

Removed useless computations at the end of procedure "leapfrog".

Merged common block "matrfil" into module "parafilt".
1 SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
2 !
3 ! Auteurs: P.Le Van, F.Hourdin, F.Forget
4 !
5 ! ***********************************************************
6 ! Shema d'advection " pseudo amont " .
7 ! *************************************************************
8 ! q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
9 ! dq sont des arguments de sortie pour le s-pg ....
10 !
11 !
12 ! ----------------------------------------------------------------
13 use dimens_m
14 use paramet_m
15 use comconst
16 use comvert
17 use logic
18 use comgeom
19 IMPLICIT NONE
20 !
21 !
22 !
23 ! Arguments:
24 ! ----------
25 REAL masse(ip1jmp1,llm),pente_max
26 REAL masse_adv_v( ip1jm,llm)
27 REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
28 !
29 ! Local
30 ! ---------
31 !
32 INTEGER i,ij,l
33 !
34 REAL airej2,airejjm,airescb(iim),airesch(iim)
35 REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
36 REAL adyqv(ip1jm),dyqmax(ip1jmp1)
37 REAL qbyv(ip1jm,llm)
38
39 REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
40 ! REAL newq,oldmasse
41 Logical extremum,first,testcpu
42 REAL temps0,temps1,temps2,temps3,temps4,temps5,second
43 SAVE temps0,temps1,temps2,temps3,temps4,temps5
44 SAVE first,testcpu
45
46 REAL convpn,convps,convmpn,convmps
47 real massepn,masseps,qpn,qps
48 REAL sinlon(iip1),sinlondlon(iip1)
49 REAL coslon(iip1),coslondlon(iip1)
50 SAVE sinlon,coslon,sinlondlon,coslondlon
51 SAVE airej2,airejjm
52 !
53 !
54 REAL SSUM
55
56 DATA first,testcpu/.true.,.false./
57 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
58
59 IF(first) THEN
60 PRINT*,'Shema Amont nouveau appele dans Vanleer '
61 first=.false.
62 do i=2,iip1
63 coslon(i)=cos(rlonv(i))
64 sinlon(i)=sin(rlonv(i))
65 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
66 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
67 ENDDO
68 coslon(1)=coslon(iip1)
69 coslondlon(1)=coslondlon(iip1)
70 sinlon(1)=sinlon(iip1)
71 sinlondlon(1)=sinlondlon(iip1)
72 airej2 = SSUM( iim, aire(iip2), 1 )
73 airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
74 ENDIF
75
76 !
77 DO l = 1, llm
78 !
79 ! --------------------------------
80 ! CALCUL EN LATITUDE
81 ! --------------------------------
82
83 ! On commence par calculer la valeur du traceur moyenne sur le
84 ! premier cercle
85 ! de latitude autour du pole (qpns pour le pole nord et qpsn pour
86 ! le pole nord) qui sera utilisee pour evaluer les pentes au pole.
87
88 DO i = 1, iim
89 airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
90 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
91 ENDDO
92 qpns = SSUM( iim, airescb ,1 ) / airej2
93 qpsn = SSUM( iim, airesch ,1 ) / airejjm
94
95 ! calcul des pentes aux points v
96
97 DO ij=1,ip1jm
98 dyqv(ij)=q(ij,l)-q(ij+iip1,l)
99 adyqv(ij)=abs(dyqv(ij))
100 ENDDO
101
102 ! calcul des pentes aux points scalaires
103
104 DO ij=iip2,ip1jm
105 dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
106 dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
107 dyqmax(ij)=pente_max*dyqmax(ij)
108 ENDDO
109
110 ! calcul des pentes aux poles
111
112 DO ij=1,iip1
113 dyq(ij,l)=qpns-q(ij+iip1,l)
114 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
115 ENDDO
116
117 ! filtrage de la derivee
118 dyn1=0.
119 dys1=0.
120 dyn2=0.
121 dys2=0.
122 DO ij=1,iim
123 dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
124 dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
125 dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
126 dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
127 ENDDO
128 DO ij=1,iip1
129 dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
130 dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
131 ENDDO
132
133 ! calcul des pentes limites aux poles
134
135 goto 8888
136 fn=1.
137 fs=1.
138 DO ij=1,iim
139 IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
140 fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
141 ENDIF
142 IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
143 fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
144 ENDIF
145 ENDDO
146 DO ij=1,iip1
147 dyq(ij,l)=fn*dyq(ij,l)
148 dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
149 ENDDO
150 8888 continue
151 DO ij=1,iip1
152 dyq(ij,l)=0.
153 dyq(ip1jm+ij,l)=0.
154 ENDDO
155
156 ! calcul des pentes limitees
157
158 DO ij=iip2,ip1jm
159 IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
160 dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
161 ELSE
162 dyq(ij,l)=0.
163 ENDIF
164 ENDDO
165
166 ENDDO
167
168 DO l=1,llm
169 DO ij=1,ip1jm
170 IF(masse_adv_v(ij,l).gt.0) THEN
171 qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)* &
172 & 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
173 ELSE
174 qbyv(ij,l)=q(ij,l)-dyq(ij,l)* &
175 & 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
176 ENDIF
177 qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
178 ENDDO
179 ENDDO
180
181
182 DO l=1,llm
183 DO ij=iip2,ip1jm
184 newmasse=masse(ij,l) &
185 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
186 q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l)) &
187 & /newmasse
188 masse(ij,l)=newmasse
189 ENDDO
190 !.-. ancienne version
191 ! convpn=SSUM(iim,qbyv(1,l),1)/apoln
192 ! convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
193
194 convpn=SSUM(iim,qbyv(1,l),1)
195 convmpn=ssum(iim,masse_adv_v(1,l),1)
196 massepn=ssum(iim,masse(1,l),1)
197 qpn=0.
198 do ij=1,iim
199 qpn=qpn+masse(ij,l)*q(ij,l)
200 enddo
201 qpn=(qpn+convpn)/(massepn+convmpn)
202 do ij=1,iip1
203 q(ij,l)=qpn
204 enddo
205
206 ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
207 ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
208
209 convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
210 convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
211 masseps=ssum(iim, masse(ip1jm+1,l),1)
212 qps=0.
213 do ij = ip1jm+1,ip1jmp1-1
214 qps=qps+masse(ij,l)*q(ij,l)
215 enddo
216 qps=(qps+convps)/(masseps+convmps)
217 do ij=ip1jm+1,ip1jmp1
218 q(ij,l)=qps
219 enddo
220
221 !.-. fin ancienne version
222
223 !._. nouvelle version
224 ! convpn=SSUM(iim,qbyv(1,l),1)
225 ! convmpn=ssum(iim,masse_adv_v(1,l),1)
226 ! oldmasse=ssum(iim,masse(1,l),1)
227 ! newmasse=oldmasse+convmpn
228 ! newq=(q(1,l)*oldmasse+convpn)/newmasse
229 ! newmasse=newmasse/apoln
230 ! DO ij = 1,iip1
231 ! q(ij,l)=newq
232 ! masse(ij,l)=newmasse*aire(ij)
233 ! ENDDO
234 ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
235 ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
236 ! oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
237 ! newmasse=oldmasse+convmps
238 ! newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
239 ! newmasse=newmasse/apols
240 ! DO ij = ip1jm+1,ip1jmp1
241 ! q(ij,l)=newq
242 ! masse(ij,l)=newmasse*aire(ij)
243 ! ENDDO
244 !._. fin nouvelle version
245 ENDDO
246
247 RETURN
248 END
  ViewVC Help
Powered by ViewVC 1.1.21