source: LMDZ5/branches/LMDZ5-DOFOCO/libf/dyn3dmem/advect_new_loc.F @ 4785

Last change on this file since 4785 was 1632, checked in by Laurent Fairhead, 12 years ago

Import initial du répertoire dyn3dmem

Attention! ceci n'est qu'une version préliminaire du code "basse mémoire":
le code contenu dans ce répertoire est basé sur la r1320 et a donc besoin
d'être mis à jour par rapport à la dynamique parallèle d'aujourd'hui.
Ce code est toutefois mis à disposition pour circonvenir à des problèmes
de mémoire que certaines configurations du modèle pourraient rencontrer.
Dans l'état, il compile et tourne sur vargas et au CCRT

Initial import of dyn3dmem

Warning! this is just a preliminary version of the memory light code:
it is based on r1320 of the code and thus needs to be updated before
it can replace the present dyn3dpar code. It is nevertheless put at your
disposal to circumvent some memory problems some LMDZ configurations may
encounter. In its present state, it will compile and run on vargas and CCRT

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
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.