source: LMDZ4/trunk/libf/dyn3dpar/advect_new_p.F @ 5440

Last change on this file since 5440 was 774, checked in by Laurent Fairhead, 18 years ago

Suite du merge entre la version et la HEAD: quelques modifications de
Yann sur le
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.0 KB
RevLine 
[774]1!
2! $Header$
3!
[763]4      SUBROUTINE advect_new_p(ucov,vcov,teta,w,massebx,masseby,
5     &                        du,dv,dteta)
6      USE parallel
7      USE write_field_p
8      IMPLICIT NONE
9c=======================================================================
10c
11c   Auteurs:  P. Le Van , Fr. Hourdin  .
12c   -------
13c
14c   Objet:
15c   ------
16c
17c   *************************************************************
18c   .... calcul des termes d'advection vertic.pour u,v,teta,q ...
19c   *************************************************************
20c        ces termes sont ajoutes a du,dv,dteta et dq .
21c  Modif F.Forget 03/94 : on retire q de advect
22c
23c=======================================================================
24c-----------------------------------------------------------------------
25c   Declarations:
26c   -------------
27
28#include "dimensions.h"
29#include "paramet.h"
30#include "comconst.h"
31#include "comvert.h"
32#include "comgeom.h"
33#include "logic.h"
34#include "ener.h"
35
36c   Arguments:
37c   ----------
38
39      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
40      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm),w(ip1jmp1,llm)
41      REAL dv(ip1jm,llm),du(ip1jmp1,llm),dteta(ip1jmp1,llm)
42      REAL,SAVE :: dv1(ip1jm,llm),du1(ip1jmp1,llm),dteta1(ip1jmp1,llm)
43      REAL,SAVE :: dv2(ip1jm,llm),du2(ip1jmp1,llm),dteta2(ip1jmp1,llm)
44c   Local:
45c   ------
46
47      REAL,SAVE :: uav(ip1jmp1,llm),vav(ip1jm,llm)
48      REAL wsur2(ip1jmp1)
49      REAL unsaire2(ip1jmp1), ge(ip1jmp1)
50      REAL deuxjour, ww, gt, uu, vv
51
52      INTEGER  ij,l,ijb,ije
53
54      EXTERNAL  SSUM
55      REAL      SSUM
56
57c-----------------------------------------------------------------------
58c   2. Calculs preliminaires:
59c   -------------------------
60
61      IF (conser)  THEN
62         deuxjour = 2. * daysec
63
64         DO   1  ij   = 1, ip1jmp1
65         unsaire2(ij) = unsaire(ij) * unsaire(ij)
66   1     CONTINUE
67      END IF
68
69
70c------------------  -yy ----------------------------------------------
71c   .  Calcul de     u
72
73c$OMP MASTER
74      ijb=ij_begin
75      ije=ij_end
76      if (pole_nord) ijb=ijb+iip1
77      if (pole_sud)  ije=ije-iip1
78
79      DO ij=ijb,ije
80        du2(ij,1)=0.
81      ENDDO
82     
83      ijb=ij_begin
84      ije=ij_end
85      if (pole_sud)  ije=ij_end-iip1
86     
87      DO ij=ijb,ije
88        dv2(ij,1)=0.
89      ENDDO
90     
91      ijb=ij_begin
92      ije=ij_end
93
94      DO ij=ijb,ije
95        dteta2(ij,1)=0.
96      ENDDO
97c$OMP END MASTER
98
99 
100c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
101      DO  l=1,llm
102         
103         ijb=ij_begin
104         ije=ij_end
105         if (pole_nord) ijb=ijb+iip1
106         if (pole_sud)  ije=ije-iip1
107         
108c         DO    ij     = iip2, ip1jmp1
109c            uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) )
110c         ENDDO
111
112c         DO    ij     = iip2, ip1jm
113c            uav(ij,l) = uav(ij,l) + uav(ij+iip1,l)
114c         ENDDO
115         
116         DO    ij     = ijb, ije
117                 
118           uav(ij,l)=0.25*(ucov(ij,l)+ucov(ij-iip1,l))
119     .               +0.25*(ucov(ij+iip1,l)+ucov(ij,l))
120         ENDDO
121         
122         if (pole_nord) then
123           DO      ij         = 1, iip1
124              uav(ij      ,l) = 0.
125           ENDDO
126         endif
127         
128         if (pole_sud) then
129           DO      ij         = 1, iip1
130              uav(ip1jm+ij,l) = 0.
131           ENDDO
132         endif
133         
134      ENDDO
135c$OMP END DO     
136c      call write_field3d_p('uav',reshape(uav,(/iip1,jjp1,llm/)))
137     
138c------------------  -xx ----------------------------------------------
139c   .  Calcul de     v
140     
141      ijb=ij_begin
142      ije=ij_end
143      if (pole_sud)  ije=ij_end-iip1
144
145c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
146      DO  l=1,llm
147         
148         DO    ij   = ijb+1, ije
149           vav(ij,l) = 0.25 * ( vcov(ij,l) + vcov(ij-1,l) )
150         ENDDO
151         
152         DO    ij   = ijb,ije,iip1
153          vav(ij,l) = vav(ij+iim,l)
154         ENDDO
155         
156         
157         DO    ij   = ijb, ije-1
158          vav(ij,l) = vav(ij,l) + vav(ij+1,l)
159         ENDDO
160         
161         DO    ij       = ijb, ije, iip1
162          vav(ij+iim,l) = vav(ij,l)
163         ENDDO
164         
165      ENDDO
166c$OMP END DO
167c       call write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/)))
168
169c-----------------------------------------------------------------------
170c$OMP BARRIER
171
172c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
173      DO 20 l = 1, llmm1
174
175
176c       ......   calcul de  - w/2.    au niveau  l+1   .......
177      ijb=ij_begin
178      ije=ij_end+iip1
179      if (pole_sud)  ije=ij_end
180     
181      DO 5   ij   = ijb, ije
182      wsur2( ij ) = - 0.5 * w( ij,l+1 )
183   5  CONTINUE
184
185
186c     .....................     calcul pour  du     ..................
187     
188      ijb=ij_begin
189      ije=ij_end
190      if (pole_nord) ijb=ijb+iip1
191      if (pole_sud)  ije=ije-iip1
192         
193      DO 6 ij = ijb ,ije-1
194      ww        = wsur2 (  ij  )     + wsur2( ij+1 )
195      uu        = 0.5 * ( ucov(ij,l) + ucov(ij,l+1) )
196      du1(ij,l)  =  ww * ( uu - uav(ij, l ) )/massebx(ij, l )
197      du2(ij,l+1)=  ww * ( uu - uav(ij,l+1) )/massebx(ij,l+1)
198   6  CONTINUE
199
200c     .................    calcul pour   dv      .....................
201      ijb=ij_begin
202      ije=ij_end
203      if (pole_sud)  ije=ij_end-iip1
204     
205      DO 8 ij = ijb, ije
206      ww        = wsur2( ij+iip1 )   + wsur2( ij )
207      vv        = 0.5 * ( vcov(ij,l) + vcov(ij,l+1) )
208      dv1(ij,l)  =  ww * (vv - vav(ij, l ) )/masseby(ij, l )
209      dv2(ij,l+1)=  ww * (vv - vav(ij,l+1) )/masseby(ij,l+1)
210   8  CONTINUE
211
212c
213
214c     ............................................................
215c     ...............    calcul pour   dh      ...................
216c     ............................................................
217
218c                       ---z
219c       calcul de  - d( teta  * w )      qu'on ajoute a   dh
220c                   ...............
221        ijb=ij_begin
222        ije=ij_end
223       
224        DO 15 ij = ijb, ije
225         ww            = wsur2(ij) * (teta(ij,l) + teta(ij,l+1) )
226         dteta1(ij, l ) =   ww
227         dteta2(ij,l+1) =   ww
228  15    CONTINUE
229
230c ym ---> conser a voir plus tard
231
232c      IF( conser)  THEN
233c       
234c        DO 17 ij = 1,ip1jmp1
235c        ge(ij)   = wsur2(ij) * wsur2(ij) * unsaire2(ij)
236c  17    CONTINUE
237c        gt       = SSUM( ip1jmp1,ge,1 )
238c        gtot(l)  = deuxjour * SQRT( gt/ip1jmp1 )
239c      END IF
240
241  20  CONTINUE
242c$OMP END DO
243
244      ijb=ij_begin
245      ije=ij_end
246      if (pole_nord) ijb=ijb+iip1
247      if (pole_sud)  ije=ije-iip1
248     
249c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
250      DO l=1,llm
251        DO ij=ijb,ije-1
252          du(ij,l)=du(ij,l)+du2(ij,l)-du1(ij,l)
253        ENDDO
254
255        DO   ij   = ijb+iip1-1, ije, iip1
256         du( ij, l  ) = du( ij -iim, l  )
257        ENDDO
258      ENDDO
259c$OMP END DO NOWAIT
260     
261      ijb=ij_begin
262      ije=ij_end
263      if (pole_sud)  ije=ij_end-iip1
264
265c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
266      DO l=1,llm
267        DO ij=ijb,ije
268          dv(ij,l)=dv(ij,l)+dv2(ij,l)-dv1(ij,l)
269        ENDDO
270      ENDDO
271c$OMP END DO NOWAIT     
272      ijb=ij_begin
273      ije=ij_end
274
275c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
276      DO l=1,llm
277        DO ij=ijb,ije
278          dteta(ij,l)=dteta(ij,l)+dteta2(ij,l)-dteta1(ij,l)
279        ENDDO
280      ENDDO
281c$OMP END DO NOWAIT     
282
283      RETURN
284      END
Note: See TracBrowser for help on using the repository browser.