source: LMDZ5/branches/LMDZ5_SPLA/libf/dyn3dmem/advect_new_loc.F

Last change on this file was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

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