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

Annotation of /trunk/dyn3d/limy.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 39 - (hide annotations)
Tue Jan 25 15:11:05 2011 UTC (13 years, 4 months ago) by guez
Original Path: trunk/libf/dyn3d/limy.f
File size: 4726 byte(s)
"pi" comes from "nr_util". Removed subroutine "initialize" in module
"comconst".

Copied the content of "fxy_sin.h" into "fxysinus", instead of getting
it from an "include" line. Removed file "fxy_sin.h".

"ps" has rank 2 in "gcm" and "dynetat0".

Assumed-shape for argument "q" of "integrd".

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 guez 39 USE nr_util, ONLY : pi
23 guez 3 IMPLICIT NONE
24     c
25     c
26     c
27     c Arguments:
28     c ----------
29     real pente_max
30     real s0(ip1jmp1,llm),sy(ip1jmp1,llm),sm(ip1jmp1,llm)
31     c
32     c Local
33     c ---------
34     c
35     INTEGER i,ij,l
36     c
37     REAL q(ip1jmp1,llm)
38     REAL airej2,airejjm,airescb(iim),airesch(iim)
39     real sigv,dyq(ip1jmp1),dyqv(ip1jm)
40     real adyqv(ip1jm),dyqmax(ip1jmp1)
41     REAL qbyv(ip1jm,llm)
42    
43     REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2
44     Logical extremum,first
45     save first
46    
47     real convpn,convps,convmpn,convmps
48     real sinlon(iip1),sinlondlon(iip1)
49     real coslon(iip1),coslondlon(iip1)
50     save sinlon,coslon,sinlondlon,coslondlon
51     c
52     c
53     REAL SSUM
54     integer ismax,ismin
55     EXTERNAL SSUM, convflu,ismin,ismax
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 cas ou on a un extremum au pole
119    
120     c if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
121     c & apn=0.
122     c if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
123     c & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
124     c & aps=0.
125    
126     c limitation des pentes aux poles
127     c do ij=1,iip1
128     c dyq(ij)=apn*dyq(ij)
129     c dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
130     c enddo
131    
132     c test
133     c do ij=1,iip1
134     c dyq(iip1+ij)=0.
135     c dyq(ip1jm+ij-iip1)=0.
136     c enddo
137     c do ij=1,ip1jmp1
138     c dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
139     c enddo
140    
141     if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
142     & then
143     do ij=1,iip1
144     dyqmax(ij)=0.
145     enddo
146     else
147     do ij=1,iip1
148     dyqmax(ij)=pente_max*abs(dyqv(ij))
149     enddo
150     endif
151    
152     if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
153     & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
154     &then
155     do ij=ip1jm+1,ip1jmp1
156     dyqmax(ij)=0.
157     enddo
158     else
159     do ij=ip1jm+1,ip1jmp1
160     dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
161     enddo
162     endif
163    
164     c calcul des pentes limitees
165    
166     do ij=1,ip1jmp1
167     if(dyqv(ij)*dyqv(ij-iip1).gt.0.) then
168     dyq(ij)=sign(min(abs(dyq(ij)),dyqmax(ij)),dyq(ij))
169     else
170     dyq(ij)=0.
171     endif
172     enddo
173    
174     DO ij=1,ip1jmp1
175     sy(ij,l) = dyq(ij) * sm ( ij,l )
176     ENDDO
177    
178     enddo ! fin de la boucle sur les couches verticales
179    
180     RETURN
181     END

  ViewVC Help
Powered by ViewVC 1.1.21