source: LMDZ4/branches/V3_test/libf/dyn3dpar/advect_new_p.F @ 727

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

Versions parallèlisées des routines YM
LF

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