My Project
 All Classes Files Functions Variables Macros
limy.F
Go to the documentation of this file.
1 c
2 c $Id: limy.F 1520 2011-05-23 11:37:09Z emillour $
3 c
4  SUBROUTINE limy(s0,sy,sm,pente_max)
5 c
6 c Auteurs: P.Le Van, F.Hourdin, F.Forget
7 c
8 c ********************************************************************
9 c Shema d'advection " pseudo amont " .
10 c ********************************************************************
11 c q,w sont des arguments d'entree pour le s-pg ....
12 c dq sont des arguments de sortie pour le s-pg ....
13 c
14 c
15 c --------------------------------------------------------------------
16  IMPLICIT NONE
17 c
18 #include "dimensions.h"
19 #include "paramet.h"
20 #include "logic.h"
21 #include "comvert.h"
22 #include "comconst.h"
23 #include "comgeom.h"
24 c
25 c
26 c Arguments:
27 c ----------
28  real pente_max
29  real s0(ip1jmp1,llm),sy(ip1jmp1,llm),sm(ip1jmp1,llm)
30 c
31 c Local
32 c ---------
33 c
34  INTEGER i,ij,l
35 c
36  REAL q(ip1jmp1,llm)
37  REAL airej2,airejjm,airescb(iim),airesch(iim)
38  real sigv,dyq(ip1jmp1),dyqv(ip1jm)
39  real adyqv(ip1jm),dyqmax(ip1jmp1)
40  REAL qbyv(ip1jm,llm)
41 
42  REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2
43  Logical extremum,first
44  save first
45 
46  real convpn,convps,convmpn,convmps
47  real sinlon(iip1),sinlondlon(iip1)
48  real coslon(iip1),coslondlon(iip1)
49  save sinlon,coslon,sinlondlon,coslondlon
50 c
51 c
52  REAL ssum
53  integer ismax,ismin
54  EXTERNAL ssum, convflu,ismin,ismax
55  EXTERNAL filtreg
56 
57  data first/.true./
58 
59  if(first) then
60  print*,'SCHEMA AMONT NOUVEAU'
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  endif
73 
74 c
75 
76  do l = 1, llm
77 c
78  DO ij=1,ip1jmp1
79  q(ij,l) = s0(ij,l) / sm( ij,l )
80  dyq(ij) = sy(ij,l) / sm( ij,l )
81  ENDDO
82 c
83 c --------------------------------
84 c CALCUL EN LATITUDE
85 c --------------------------------
86 
87 c On commence par calculer la valeur du traceur moyenne sur le premier cercle
88 c de latitude autour du pole (qpns pour le pole nord et qpsn pour
89 c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
90 
91  airej2 = ssum( iim, aire(iip2), 1 )
92  airejjm= ssum( iim, aire(ip1jm -iim), 1 )
93  DO i = 1, iim
94  airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
95  airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
96  ENDDO
97  qpns = ssum( iim, airescb ,1 ) / airej2
98  qpsn = ssum( iim, airesch ,1 ) / airejjm
99 
100 c calcul des pentes aux points v
101 
102  do ij=1,ip1jm
103  dyqv(ij)=q(ij,l)-q(ij+iip1,l)
104  adyqv(ij)=abs(dyqv(ij))
105  ENDDO
106 
107 c calcul des pentes aux points scalaires
108 
109  do ij=iip2,ip1jm
110  dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
111  dyqmax(ij)=pente_max*dyqmax(ij)
112  enddo
113 
114 c calcul des pentes aux poles
115 
116 c calcul des pentes limites aux poles
117 
118 c print*,dyqv(iip1+1)
119 c appn=abs(dyq(1)/dyqv(iip1+1))
120 c print*,dyq(ip1jm+1)
121 c print*,dyqv(ip1jm-iip1+1)
122 c apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
123 c do ij=2,iim
124 c appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
125 c apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
126 c enddo
127 c appn=min(pente_max/appn,1.)
128 c apps=min(pente_max/apps,1.)
129 
130 
131 c cas ou on a un extremum au pole
132 
133 c if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
134 c & appn=0.
135 c if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
136 c & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
137 c & apps=0.
138 
139 c limitation des pentes aux poles
140 c do ij=1,iip1
141 c dyq(ij)=appn*dyq(ij)
142 c dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
143 c enddo
144 
145 c test
146 c do ij=1,iip1
147 c dyq(iip1+ij)=0.
148 c dyq(ip1jm+ij-iip1)=0.
149 c enddo
150 c do ij=1,ip1jmp1
151 c dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
152 c enddo
153 
154  if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
155  & then
156  do ij=1,iip1
157  dyqmax(ij)=0.
158  enddo
159  else
160  do ij=1,iip1
161  dyqmax(ij)=pente_max*abs(dyqv(ij))
162  enddo
163  endif
164 
165  if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
166  & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
167  &then
168  do ij=ip1jm+1,ip1jmp1
169  dyqmax(ij)=0.
170  enddo
171  else
172  do ij=ip1jm+1,ip1jmp1
173  dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
174  enddo
175  endif
176 
177 c calcul des pentes limitees
178 
179  do ij=1,ip1jmp1
180  if(dyqv(ij)*dyqv(ij-iip1).gt.0.) then
181  dyq(ij)=sign(min(abs(dyq(ij)),dyqmax(ij)),dyq(ij))
182  else
183  dyq(ij)=0.
184  endif
185  enddo
186 
187  DO ij=1,ip1jmp1
188  sy(ij,l) = dyq(ij) * sm( ij,l )
189  ENDDO
190 
191  enddo ! fin de la boucle sur les couches verticales
192 
193  RETURN
194  END