source: codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90 @ 362

Last change on this file since 362 was 362, checked in by dubos, 9 years ago

Introduced DEC convention for velocity - HEVI only - cleanup to follow
File size: 19.9 KB
Line 
1MODULE caldyn_kernels_hevi_mod
2  USE icosa
3  USE transfert_mod
4  USE caldyn_kernels_base_mod
5  IMPLICIT NONE
6  PRIVATE
7
8  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_fast, compute_caldyn_slow
9
10CONTAINS
11
12  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
13    USE icosa
14    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop
15    USE exner_mod
16    USE trace
17    USE omp_para
18    IMPLICIT NONE
19    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
20    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
21    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
22    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
23    INTEGER :: ij,l
24    REAL(rstd) :: m
25    CALL trace_start("compute_theta") 
26
27    IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
28       DO l = ll_begin,ll_end
29          !DIR$ SIMD
30          DO ij=ij_begin_ext,ij_end_ext
31             IF(DEC) THEN  ! ps is actually Ms
32                m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l)
33             ELSE
34                m = mass_dak(l)+ps(ij)*mass_dbk(l)
35             END IF
36             rhodz(ij,l) = m/g
37             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
38          ENDDO
39       ENDDO
40    ELSE ! Compute only theta
41       DO l = ll_begin,ll_end
42          !DIR$ SIMD
43          DO ij=ij_begin_ext,ij_end_ext
44             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
45          ENDDO
46       ENDDO
47    END IF
48
49    CALL trace_end("compute_theta")
50  END SUBROUTINE compute_theta
51
52  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
53    USE icosa
54    USE exner_mod
55    USE trace
56    USE omp_para
57    IMPLICIT NONE
58    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
59    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
60    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
61    REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
62
63    INTEGER :: ij,l
64    REAL(rstd) :: etav,hv,radius_m2
65
66    CALL trace_start("compute_pvort_only") 
67!!! Compute shallow-water potential vorticity
68    radius_m2=radius**(-2)
69    DO l = ll_begin,ll_end
70       !DIR$ SIMD
71       DO ij=ij_begin_ext,ij_end_ext
72          IF(DEC) THEN
73             etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)         &
74                                   + ne_left * u(ij+t_rup+u_left,l)  &
75                                   - ne_lup  * u(ij+u_lup,l) )                               
76
77             hv =   Riv2(ij,vup)          * rhodz(ij,l)           &
78                  + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
79                  + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
80             qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
81
82             etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          &
83                                      + ne_right * u(ij+t_ldown+u_right,l)  &
84                                      - ne_rdown * u(ij+u_rdown,l) )
85             hv =   Riv2(ij,vdown)        * rhodz(ij,l)          &
86                  + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
87                  + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
88             qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
89          ELSE
90             etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)        * de(ij+u_rup)         &
91                                   + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
92                                   - ne_lup  * u(ij+u_lup,l)        * de(ij+u_lup) )                               
93
94             hv =   Riv2(ij,vup)          * rhodz(ij,l)           &
95                  + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
96                  + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
97             qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
98
99             etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
100                                      + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
101                                      - ne_rdown * u(ij+u_rdown,l)          * de(ij+u_rdown) )
102             hv =   Riv2(ij,vdown)        * rhodz(ij,l)          &
103                  + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
104                  + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
105             qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
106          END IF
107       ENDDO
108
109       !DIR$ SIMD
110       DO ij=ij_begin,ij_end
111          qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l)) 
112          qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l)) 
113          qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l)) 
114       END DO
115
116    ENDDO
117
118    CALL trace_end("compute_pvort_only")
119
120  END SUBROUTINE compute_pvort_only
121
122  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot, du)
123    USE icosa
124    USE disvert_mod
125    USE exner_mod
126    USE trace
127    USE omp_para
128    IMPLICIT NONE
129    REAL(rstd), INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs
130    REAL(rstd),INTENT(INOUT)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
131    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
132    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
133    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
134    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
135    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
136    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
137
138    INTEGER :: i,j,ij,l
139    REAL(rstd) :: due_right, due_lup, due_ldown
140
141    CALL trace_start("compute_caldyn_fast")
142   
143    ! Compute bernouilli term
144    IF(boussinesq) THEN
145       DO l=ll_begin,ll_end
146          !DIR$ SIMD
147          DO ij=ij_begin,ij_end         
148             berni(ij,l) = pk(ij,l)
149             ! from now on pk contains the vertically-averaged geopotential
150             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
151          ENDDO
152       ENDDO
153
154    ELSE ! compressible
155
156       DO l=ll_begin,ll_end
157          !DIR$ SIMD
158          DO ij=ij_begin,ij_end         
159             berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
160          ENDDO
161       ENDDO
162
163    END IF ! Boussinesq/compressible
164
165!!! u:=u+tau*du, du = gradients of Bernoulli and Exner functions
166    DO l=ll_begin,ll_end
167       !DIR$ SIMD
168       DO ij=ij_begin,ij_end         
169          due_right =  0.5*(theta(ij,l)+theta(ij+t_right,l))                  &
170                          *(ne_right*pk(ij,l)   +ne_left*pk(ij+t_right,l))    &
171                         +  ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l)
172          due_lup = 0.5*(theta(ij,l)+theta(ij+t_lup,l))                       &
173                       *(ne_lup*pk(ij,l)   +ne_rdown*pk(ij+t_lup,l))          &
174                      +  ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l)
175          due_ldown = 0.5*(theta(ij,l)+theta(ij+t_ldown,l))                   &
176                         *(ne_ldown*pk(ij,l)   +ne_rup*pk(ij+t_ldown,l))      &
177                        +  ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l)
178          IF(.NOT.DEC) THEN
179             due_right = due_right/de(ij+u_right)
180             due_lup   = due_lup/de(ij+u_lup)
181             due_ldown = due_ldown/de(ij+u_ldown)
182          END IF
183          du(ij+u_right,l) = due_right
184          du(ij+u_lup,l)   = due_lup
185          du(ij+u_ldown,l) = due_ldown
186          u(ij+u_right,l) = u(ij+u_right,l) + tau*due_right
187          u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*due_lup
188          u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*due_ldown
189       ENDDO
190    ENDDO
191
192    CALL trace_end("compute_caldyn_fast")
193
194  END SUBROUTINE compute_caldyn_fast
195
196  SUBROUTINE compute_caldyn_slow(u,rhodz,qu,theta, hflux,convm, dtheta_rhodz, du)
197    USE icosa
198    USE disvert_mod
199    USE exner_mod
200    USE trace
201    USE omp_para
202    IMPLICIT NONE
203    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
204    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
205    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
206    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
207
208    REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s
209    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
210    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
211    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
212
213    REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi)
214    REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
215    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
216    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
217
218    INTEGER :: ij,l
219    REAL(rstd) :: uu_right, uu_lup, uu_ldown
220
221    CALL trace_start("compute_caldyn_slow")
222
223    DO l = ll_begin, ll_end
224!!!  Compute mass and theta fluxes
225       IF (caldyn_conserv==energy) CALL test_message(req_qu) 
226       !DIR$ SIMD
227       DO ij=ij_begin_ext,ij_end_ext
228          uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
229          uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
230          uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
231          IF(DEC) THEN
232             uu_right= uu_right*le(ij+u_right)/de(ij+u_right)
233             uu_lup  = uu_lup  *le(ij+u_lup)/de(ij+u_lup)
234             uu_ldown= uu_ldown*le(ij+u_ldown)/de(ij+u_ldown)
235          ELSE
236             uu_right= uu_right*le(ij+u_right)
237             uu_lup  = uu_lup  *le(ij+u_lup)
238             uu_ldown= uu_ldown*le(ij+u_ldown)
239          END IF
240          hflux(ij+u_right,l)=uu_right
241          hflux(ij+u_lup,l)  =uu_lup
242          hflux(ij+u_ldown,l)=uu_ldown
243          Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*uu_right
244          Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*uu_lup
245          Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*uu_ldown
246       ENDDO
247
248!!! compute horizontal divergence of fluxes
249       !DIR$ SIMD
250       DO ij=ij_begin,ij_end         
251          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
252          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
253               ne_rup*hflux(ij+u_rup,l)       +  & 
254               ne_lup*hflux(ij+u_lup,l)       +  & 
255               ne_left*hflux(ij+u_left,l)     +  & 
256               ne_ldown*hflux(ij+u_ldown,l)   +  & 
257               ne_rdown*hflux(ij+u_rdown,l))   
258
259          ! signe ? attention d (rho theta dz)
260          ! dtheta_rhodz =  -div(flux.theta)
261          dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
262               ne_rup*Ftheta(ij+u_rup,l)        +  & 
263               ne_lup*Ftheta(ij+u_lup,l)        +  & 
264               ne_left*Ftheta(ij+u_left,l)      +  & 
265               ne_ldown*Ftheta(ij+u_ldown,l)    +  & 
266               ne_rdown*Ftheta(ij+u_rdown,l))
267       ENDDO
268
269    END DO
270
271!!! Compute potential vorticity (Coriolis) contribution to du
272
273    SELECT CASE(caldyn_conserv)
274    CASE(energy) ! energy-conserving TRiSK
275
276       CALL wait_message(req_qu)
277
278       DO l=ll_begin,ll_end
279          !DIR$ SIMD
280          DO ij=ij_begin,ij_end         
281             uu_right = &
282                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
283                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
284                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
285                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
286                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
287                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
288                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
289                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
290                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
291                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
292             uu_lup = &
293                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
294                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
295                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
296                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
297                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
298                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
299                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
300                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
301                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
302                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
303             uu_ldown = &
304                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
305                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
306                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
307                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
308                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
309                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
310                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
311                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
312                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
313                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
314             IF(DEC) THEN
315                du(ij+u_right,l) = .5*uu_right
316                du(ij+u_lup,l)   = .5*uu_lup
317                du(ij+u_ldown,l)   = .5*uu_ldown
318             ELSE
319                du(ij+u_right,l) = .5*uu_right/de(ij+u_right)
320                du(ij+u_lup,l)   = .5*uu_lup  /de(ij+u_lup)
321                du(ij+u_ldown,l) = .5*uu_ldown/de(ij+u_ldown)
322             END IF
323          ENDDO
324       ENDDO
325
326    CASE(enstrophy) ! enstrophy-conserving TRiSK
327
328       DO l=ll_begin,ll_end
329          !DIR$ SIMD
330          DO ij=ij_begin,ij_end         
331             uu_right = &
332                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
333                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
334                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
335                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
336                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
337                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
338                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
339                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
340                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
341                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
342             uu_lup = &
343                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
344                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
345                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
346                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
347                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
348                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
349                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
350                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
351                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
352                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
353             uu_ldown = &
354                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
355                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
356                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
357                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
358                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
359                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
360                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
361                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
362                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
363                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
364             IF(DEC) THEN
365                du(ij+u_right,l) = qu(ij+u_right,l)*uu_right
366                du(ij+u_lup,l)   = qu(ij+u_lup,l)  *uu_lup
367                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu_ldown
368             ELSE
369                du(ij+u_right,l) = qu(ij+u_right,l)*uu_right/de(ij+u_right)
370                du(ij+u_lup,l)   = qu(ij+u_lup,l)  *uu_lup  /de(ij+u_lup)
371                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu_ldown/de(ij+u_ldown)
372             END IF
373          ENDDO
374       ENDDO
375
376    CASE DEFAULT
377       STOP
378    END SELECT
379
380    ! Compute bernouilli term = Kinetic Energy
381    le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect
382    DO l=ll_begin,ll_end
383       !DIR$ SIMD
384       DO ij=ij_begin,ij_end         
385          IF(DEC) THEN
386             berni(ij,l) = &
387                  1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
388                  le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
389                  le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
390                  le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
391                  le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
392                  le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
393          ELSE
394             berni(ij,l) = &
395                  1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
396                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
397                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
398                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
399                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
400                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
401          END IF
402       ENDDO
403    ENDDO
404
405!!! Add gradients of Bernoulli and Exner functions to du
406    DO l=ll_begin,ll_end
407       !DIR$ SIMD
408       DO ij=ij_begin,ij_end
409          IF(DEC) THEN
410             du(ij+u_right,l) = du(ij+u_right,l)                       &
411                  + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l)
412             du(ij+u_lup,l) = du(ij+u_lup,l)                           &
413                  + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l)
414             du(ij+u_ldown,l) = du(ij+u_ldown,l)                       &
415                  + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l)
416          ELSE
417             du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right)            &
418                  * ( ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
419             du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup)                 &
420                  * ( ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
421             du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown)            &
422                  * ( ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
423          END IF
424       END DO
425    ENDDO
426
427    CALL trace_end("compute_caldyn_slow")
428
429  END SUBROUTINE compute_caldyn_slow
430
431END MODULE caldyn_kernels_hevi_mod
Note: See TracBrowser for help on using the repository browser.