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

Annotation of /trunk/dyn3d/Vlsplt/vly.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 32 - (hide 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 guez 31 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