source: LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/limy.F @ 4667

Last change on this file since 4667 was 524, checked in by lmdzadmin, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
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,apn,aps,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, 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
74c
75
76      do l = 1, llm
77c
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
82c
83c   --------------------------------
84c      CALCUL EN LATITUDE
85c   --------------------------------
86
87c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
88c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
89c    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
100c   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
107c   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
114c   calcul des pentes aux poles
115
116c   calcul des pentes limites aux poles
117
118c     print*,dyqv(iip1+1)
119c     apn=abs(dyq(1)/dyqv(iip1+1))
120c     print*,dyq(ip1jm+1)
121c     print*,dyqv(ip1jm-iip1+1)
122c     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
123c     do ij=2,iim
124c        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
125c        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
126c     enddo
127c     apn=min(pente_max/apn,1.)
128c     aps=min(pente_max/aps,1.)
129
130
131c   cas ou on a un extremum au pole
132
133c     if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
134c    &   apn=0.
135c     if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
136c    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
137c    &   aps=0.
138
139c   limitation des pentes aux poles
140c     do ij=1,iip1
141c        dyq(ij)=apn*dyq(ij)
142c        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
143c     enddo
144
145c   test
146c      do ij=1,iip1
147c         dyq(iip1+ij)=0.
148c         dyq(ip1jm+ij-iip1)=0.
149c      enddo
150c      do ij=1,ip1jmp1
151c         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
152c      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
177c   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
Note: See TracBrowser for help on using the repository browser.