source: LMDZ6/trunk/libf/dyn3dmem/advect_new_loc.f90 @ 5281

Last change on this file since 5281 was 5281, checked in by abarral, 4 days ago

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