/[lmdze]/trunk/libf/dyn3d/vlyqs.f90
ViewVC logotype

Contents of /trunk/libf/dyn3d/vlyqs.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 43 - (show annotations)
Fri Apr 8 12:43:31 2011 UTC (13 years, 1 month ago) by guez
File size: 6755 byte(s) 
"start_init_phys" is now called directly by "etat0" instead of through
"start_init_dyn". "qsol_2d" is no longer a variable of module
"start_init_phys_m", it is an argument of
"start_init_phys". "start_init_dyn" now receives "tsol_2d" from
"etat0".

Split file "vlspltqs.f" into "vlspltqs.f90", "vlxqs.f90" and
""vlyqs.f90".

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