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

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

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F to *.f90
Fix gradsdef.h formatting
Remove unnecessary "RETURN" at the end of functions/subroutines

  • 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.8 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  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    DO    ij = ijb, ije - 1
161      vav(ij, l) = vav(ij, l) + vav(ij + 1, l)
162    ENDDO
163
164    DO    ij = ijb, ije, iip1
165      vav(ij + iim, l) = vav(ij, l)
166    ENDDO
167
168  ENDDO
169  !$OMP END DO
170  ! CALL write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/)))
171
172  !-----------------------------------------------------------------------
173  !$OMP BARRIER
174
175  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
176  DO l = 1, llmm1
177
178
179    ! ......   calcul de  - w/2.    au niveau  l+1   .......
180    ijb = ij_begin
181    ije = ij_end + iip1
182    if (pole_sud)  ije = ij_end
183
184    DO ij = ijb, ije
185      wsur2(ij) = - 0.5 * w(ij, l + 1)
186    END DO
187
188
189    ! .....................     calcul pour  du     ..................
190
191    ijb = ij_begin
192    ije = ij_end
193    if (pole_nord) ijb = ijb + iip1
194    if (pole_sud)  ije = ije - iip1
195
196    DO ij = ijb, ije - 1
197      ww = wsur2 (ij) + wsur2(ij + 1)
198      uu = 0.5 * (ucov(ij, l) + ucov(ij, l + 1))
199      du1(ij, l) = ww * (uu - uav(ij, l)) / massebx(ij, l)
200      du2(ij, l + 1) = ww * (uu - uav(ij, l + 1)) / massebx(ij, l + 1)
201    END DO
202
203    ! .................    calcul pour   dv      .....................
204    ijb = ij_begin
205    ije = ij_end
206    if (pole_sud)  ije = ij_end - iip1
207
208    DO ij = ijb, ije
209      ww = wsur2(ij + iip1) + wsur2(ij)
210      vv = 0.5 * (vcov(ij, l) + vcov(ij, l + 1))
211      dv1(ij, l) = ww * (vv - vav(ij, l)) / masseby(ij, l)
212      dv2(ij, l + 1) = ww * (vv - vav(ij, l + 1)) / masseby(ij, l + 1)
213    END DO
214
215    !
216
217    ! ............................................................
218    ! ...............    calcul pour   dh      ...................
219    ! ............................................................
220
221    !                   ---z
222    !   calcul de  - d( teta  * w )      qu'on ajoute a   dh
223    !               ...............
224    ijb = ij_begin
225    ije = ij_end
226
227    DO ij = ijb, ije
228      ww = wsur2(ij) * (teta(ij, l) + teta(ij, l + 1))
229      dteta1(ij, l) = ww
230      dteta2(ij, l + 1) = ww
231    END DO
232
233    ! ym ---> conser a voir plus tard
234
235    ! IF( conser)  THEN
236    !
237    !    DO 17 ij = 1,ip1jmp1
238    !    ge(ij)   = wsur2(ij) * wsur2(ij) * unsaire2(ij)
239    !  17    CONTINUE
240    !    gt       = SSUM( ip1jmp1,ge,1 )
241    !    gtot(l)  = deuxjour * SQRT( gt/ip1jmp1 )
242    !  END IF
243
244  END DO
245  !$OMP END DO
246
247  ijb = ij_begin
248  ije = ij_end
249  if (pole_nord) ijb = ijb + iip1
250  if (pole_sud)  ije = ije - iip1
251  IF (CPPKEY_DEBUGIO) THEN
252    CALL WriteField_u('du_bis', du)
253  END IF
254  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
255  DO l = 1, llm
256    DO ij = ijb, ije - 1
257      du(ij, l) = du(ij, l) + du2(ij, l) - du1(ij, l)
258    ENDDO
259
260    DO   ij = ijb + iip1 - 1, ije, iip1
261      du(ij, l) = du(ij - iim, l)
262    ENDDO
263  ENDDO
264  !$OMP END DO NOWAIT
265  IF (CPPKEY_DEBUGIO) THEN
266    CALL WriteField_u('du1', du1)
267    CALL WriteField_u('du2', du2)
268    CALL WriteField_u('du_bis', du)
269  END IF
270  ijb = ij_begin
271  ije = ij_end
272  if (pole_sud)  ije = ij_end - iip1
273
274  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
275  DO l = 1, llm
276    DO ij = ijb, ije
277      dv(ij, l) = dv(ij, l) + dv2(ij, l) - dv1(ij, l)
278    ENDDO
279  ENDDO
280  !$OMP END DO NOWAIT
281  ijb = ij_begin
282  ije = ij_end
283
284  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
285  DO l = 1, llm
286    DO ij = ijb, ije
287      dteta(ij, l) = dteta(ij, l) + dteta2(ij, l) - dteta1(ij, l)
288    ENDDO
289  ENDDO
290  !$OMP END DO NOWAIT
291
292
293END SUBROUTINE advect_new_loc
Note: See TracBrowser for help on using the repository browser.