/[lmdze]/trunk/libf/dyn3d/limy.f
ViewVC logotype

Annotation of /trunk/libf/dyn3d/limy.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (hide annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 2 months ago) by guez
File size: 5091 byte(s)
Initial import
1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/limy.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $
3     !
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     use dimens_m
17     use paramet_m
18     use comconst
19     use comvert
20     use logic
21     use comgeom
22     IMPLICIT NONE
23     c
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,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
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 apn=abs(dyq(1)/dyqv(iip1+1))
120     c print*,dyq(ip1jm+1)
121     c print*,dyqv(ip1jm-iip1+1)
122     c aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
123     c do ij=2,iim
124     c apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
125     c aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
126     c enddo
127     c apn=min(pente_max/apn,1.)
128     c aps=min(pente_max/aps,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 & apn=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 & aps=0.
138    
139     c limitation des pentes aux poles
140     c do ij=1,iip1
141     c dyq(ij)=apn*dyq(ij)
142     c dyq(ip1jm+ij)=aps*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

  ViewVC Help
Powered by ViewVC 1.1.21