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

Last change on this file since 5411 was 5159, checked in by abarral, 5 months ago

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