source: trunk/LMDZ.COMMON/libf/dyn3dpar/advect_p.F @ 1243

Last change on this file since 1243 was 1019, checked in by emillour, 11 years ago

Common dynamics; keep up with updates (seq and ) in LMDZ5 (up tio rev 1845):

  • General stuff:
  • makelmdz_fcm: add options -j # (compile using # threads) and -full, and to keep up

with Earth model, possibility to compile with various versions of orchidee

  • bld.cfg: adaptations to enable compiling using multiple threads
  • build_gcm: adaptations to enable compiling using multiple threads
  • makelmdz: keep up with Earth model: possibility to compile with various versions of orchidee + cosmetic changes + library directory name change
  • bibio:
  • wxios.F90 : Added for possible future use of XIOS library
  • filtrez:
  • mkl_dft_type.f90 & mkl_dfti.f90 : MKL (for MKL FFT) interface definitions
  • filtreg_mod : limit use of FFT to parallel mode
  • mod_filtre_fft.F90 & mod_filtre_fft_lov.F90 : swich to use parallel_lmdz
  • dyn3d:
  • abort_gcm.F : add things for xios
  • advtrac.F90 : minor change in CFL outputs
  • ce0l.F90 : indicesol.h is now module indice_sol_mod
  • comvert.h : cosmetic change on comments
  • gcm.F : add xios and use module indice_sol_mod (for INCA)
  • inigeom.F : move two computations outside loop
  • dyn3dpar:
  • parallel.F90 => parallel_lmdz.F90 : and change all the "use parallel" into "use parallel_lmdz" in all files in dyn3dpar
  • comvert.h : cosmetic change on comments
  • gcm.F : add xios and use module indice_sol_mod (for INCA)
  • leapfrog_p.F : add xios + correction for times in Newtonian case
  • ce0l.F90 : indicesol.h is now module indice_sol_mod
  • inigeom.F : move two computations outside loop

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