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

Last change on this file since 3553 was 1459, checked in by emillour, 10 years ago

Common dynamics: A couple of bug fixes

  • calfis[_p].F array boundaries must be explicitely specified when underlying arrays are of different sizes.
  • advect_new_p.F : missing initializations of intermediate variables at topmost layer.

EM
File size: 7.3 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
[1459]30      include "dimensions.h"
31      include "paramet.h"
32      include "comgeom.h"
[1]33
34c   Arguments:
35c   ----------
36
[1459]37      REAL,INTENT(IN) :: vcov(ip1jm,llm)
38      REAL,INTENT(IN) :: ucov(ip1jmp1,llm)
39      REAL,INTENT(IN) :: teta(ip1jmp1,llm)
40      REAL,INTENT(IN) :: massebx(ip1jmp1,llm)
41      REAL,INTENT(IN) :: masseby(ip1jm,llm)
42      REAL,INTENT(IN) :: w(ip1jmp1,llm)
43      REAL,INTENT(INOUT) :: dv(ip1jm,llm)
44      REAL,INTENT(INOUT) :: du(ip1jmp1,llm)
45      REAL,INTENT(INOUT) :: dteta(ip1jmp1,llm)
[1]46c   Local:
47c   ------
48
[1459]49      REAL,SAVE :: dv1(ip1jm,llm),du1(ip1jmp1,llm),dteta1(ip1jmp1,llm)
50      REAL,SAVE :: dv2(ip1jm,llm),du2(ip1jmp1,llm),dteta2(ip1jmp1,llm)
[1]51      REAL,SAVE :: uav(ip1jmp1,llm),vav(ip1jm,llm)
52      REAL wsur2(ip1jmp1)
53      REAL unsaire2(ip1jmp1), ge(ip1jmp1)
54      REAL deuxjour, ww, gt, uu, vv
55
56      INTEGER  ij,l,ijb,ije
57
58      EXTERNAL  SSUM
59      REAL      SSUM
60
61c-----------------------------------------------------------------------
62c   2. Calculs preliminaires:
63c   -------------------------
64
65      IF (conser)  THEN
66         deuxjour = 2. * daysec
67
[1459]68         DO ij   = 1, ip1jmp1
[1]69         unsaire2(ij) = unsaire(ij) * unsaire(ij)
[1459]70         ENDDO
[1]71      END IF
72
73
74c------------------  -yy ----------------------------------------------
75c   .  Calcul de     u
76
77c$OMP MASTER
78      ijb=ij_begin
79      ije=ij_end
80      if (pole_nord) ijb=ijb+iip1
81      if (pole_sud)  ije=ije-iip1
82
83      DO ij=ijb,ije
84        du2(ij,1)=0.
[1459]85        du1(ij,llm)=0.
[1]86      ENDDO
87     
88      ijb=ij_begin
89      ije=ij_end
90      if (pole_sud)  ije=ij_end-iip1
91     
92      DO ij=ijb,ije
93        dv2(ij,1)=0.
[1459]94        dv1(ij,llm)=0.
[1]95      ENDDO
96     
97      ijb=ij_begin
98      ije=ij_end
99
100      DO ij=ijb,ije
101        dteta2(ij,1)=0.
[1459]102        dteta1(ij,llm)=0.
[1]103      ENDDO
104c$OMP END MASTER
105
106 
107c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
108      DO  l=1,llm
109         
110         ijb=ij_begin
111         ije=ij_end
112         if (pole_nord) ijb=ijb+iip1
113         if (pole_sud)  ije=ije-iip1
114         
115c         DO    ij     = iip2, ip1jmp1
116c            uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) )
117c         ENDDO
118
119c         DO    ij     = iip2, ip1jm
120c            uav(ij,l) = uav(ij,l) + uav(ij+iip1,l)
121c         ENDDO
122         
123         DO    ij     = ijb, ije
124                 
125           uav(ij,l)=0.25*(ucov(ij,l)+ucov(ij-iip1,l))
126     .               +0.25*(ucov(ij+iip1,l)+ucov(ij,l))
127         ENDDO
128         
129         if (pole_nord) then
130           DO      ij         = 1, iip1
131              uav(ij      ,l) = 0.
132           ENDDO
133         endif
134         
135         if (pole_sud) then
136           DO      ij         = 1, iip1
137              uav(ip1jm+ij,l) = 0.
138           ENDDO
139         endif
[1459]140
[1]141      ENDDO
142c$OMP END DO     
143c      call write_field3d_p('uav',reshape(uav,(/iip1,jjp1,llm/)))
144     
145c------------------  -xx ----------------------------------------------
146c   .  Calcul de     v
147     
148      ijb=ij_begin
149      ije=ij_end
150      if (pole_sud)  ije=ij_end-iip1
151
152c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
153      DO  l=1,llm
154         
155         DO    ij   = ijb+1, ije
156           vav(ij,l) = 0.25 * ( vcov(ij,l) + vcov(ij-1,l) )
157         ENDDO
158         
159         DO    ij   = ijb,ije,iip1
160          vav(ij,l) = vav(ij+iim,l)
161         ENDDO
162         
163         
164         DO    ij   = ijb, ije-1
165          vav(ij,l) = vav(ij,l) + vav(ij+1,l)
166         ENDDO
167         
168         DO    ij       = ijb, ije, iip1
169          vav(ij+iim,l) = vav(ij,l)
170         ENDDO
171         
172      ENDDO
173c$OMP END DO
174c       call write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/)))
175
176c-----------------------------------------------------------------------
177c$OMP BARRIER
178
179c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1459]180      DO l = 1, llmm1
[1]181
182
183c       ......   calcul de  - w/2.    au niveau  l+1   .......
[1459]184        ijb=ij_begin
185        ije=ij_end+iip1
186        if (pole_sud)  ije=ij_end
[1]187     
[1459]188        DO ij   = ijb, ije
189          wsur2( ij ) = - 0.5 * w( ij,l+1 )
190        ENDDO
[1]191
192
193c     .....................     calcul pour  du     ..................
194     
[1459]195        ijb=ij_begin
196        ije=ij_end
197        if (pole_nord) ijb=ijb+iip1
198        if (pole_sud)  ije=ije-iip1
[1]199         
[1459]200        DO ij = ijb ,ije-1
201          ww        = wsur2 (  ij  )     + wsur2( ij+1 )
202          uu        = 0.5 * ( ucov(ij,l) + ucov(ij,l+1) )
203          du1(ij,l)  =  ww * ( uu - uav(ij, l ) )/massebx(ij, l )
204          du2(ij,l+1)=  ww * ( uu - uav(ij,l+1) )/massebx(ij,l+1)
205        ENDDO
[1]206
207c     .................    calcul pour   dv      .....................
[1459]208        ijb=ij_begin
209        ije=ij_end
210        if (pole_sud)  ije=ij_end-iip1
[1]211     
[1459]212        DO ij = ijb, ije
213          ww        = wsur2( ij+iip1 )   + wsur2( ij )
214          vv        = 0.5 * ( vcov(ij,l) + vcov(ij,l+1) )
215          dv1(ij,l)  =  ww * (vv - vav(ij, l ) )/masseby(ij, l )
216          dv2(ij,l+1)=  ww * (vv - vav(ij,l+1) )/masseby(ij,l+1)
217        ENDDO
[1]218
219c
220
221c     ............................................................
222c     ...............    calcul pour   dh      ...................
223c     ............................................................
224
225c                       ---z
226c       calcul de  - d( teta  * w )      qu'on ajoute a   dh
227c                   ...............
228        ijb=ij_begin
229        ije=ij_end
230       
[1459]231        DO ij = ijb, ije
[1]232         ww            = wsur2(ij) * (teta(ij,l) + teta(ij,l+1) )
233         dteta1(ij, l ) =   ww
234         dteta2(ij,l+1) =   ww
[1459]235        ENDDO
[1]236
237c ym ---> conser a voir plus tard
238
239c      IF( conser)  THEN
240c       
241c        DO 17 ij = 1,ip1jmp1
242c        ge(ij)   = wsur2(ij) * wsur2(ij) * unsaire2(ij)
243c  17    CONTINUE
244c        gt       = SSUM( ip1jmp1,ge,1 )
245c        gtot(l)  = deuxjour * SQRT( gt/ip1jmp1 )
246c      END IF
247
[1459]248      ENDDO ! of DO l = 1, llmm1
[1]249c$OMP END DO
250
251      ijb=ij_begin
252      ije=ij_end
253      if (pole_nord) ijb=ijb+iip1
254      if (pole_sud)  ije=ije-iip1
255     
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     
268      ijb=ij_begin
269      ije=ij_end
270      if (pole_sud)  ije=ij_end-iip1
271
272c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
273      DO l=1,llm
274        DO ij=ijb,ije
275          dv(ij,l)=dv(ij,l)+dv2(ij,l)-dv1(ij,l)
276        ENDDO
277      ENDDO
278c$OMP END DO NOWAIT     
279      ijb=ij_begin
280      ije=ij_end
281
282c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
283      DO l=1,llm
284        DO ij=ijb,ije
285          dteta(ij,l)=dteta(ij,l)+dteta2(ij,l)-dteta1(ij,l)
286        ENDDO
287      ENDDO
288c$OMP END DO NOWAIT     
289
290      END
Note: See TracBrowser for help on using the repository browser.