source: LMDZ4/branches/V3_test/libf/dyn3dpar/advect_p.F @ 5384

Last change on this file since 5384 was 630, checked in by Laurent Fairhead, 20 years ago

Import d'une version parallele de la dynamique YM
LF

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