source: trunk/LMDZ.COMMON/libf/dyn3dpar/advect_new_p.F @ 1453

Last change on this file since 1453 was 1422, checked in by milmd, 10 years ago

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