source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/advect_new_loc.f90 @ 5123

Last change on this file since 5123 was 5123, checked in by abarral, 4 months ago

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