source: trunk/libf/dyn3dpar/limy.F @ 109

Last change on this file since 109 was 109, checked in by slebonnois, 14 years ago

SLebonnois: discretisation verticale: cohabitation entre
la methode Terre et les autres.

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