source: codes/icosagcm/trunk/src/caldyn_kernels.f90 @ 361

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

ARK2.3 HEVI time stepping - tested with DCMIP4 only
File size: 42.3 KB
Line 
1MODULE caldyn_kernels_mod
2  USE icosa
3  USE transfert_mod
4  IMPLICIT NONE
5
6  INTEGER, PARAMETER :: energy=1, enstrophy=2
7  TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:)
8  REAL(rstd),SAVE,POINTER :: out_u(:,:), p(:,:), qu(:,:)
9  !$OMP THREADPRIVATE(out_u, p, qu)
10
11  ! temporary shared variable for caldyn
12  TYPE(t_field),POINTER :: f_pk(:)
13  TYPE(t_field),POINTER :: f_wwuu(:)
14  TYPE(t_field),POINTER :: f_planetvel(:)
15
16  INTEGER :: caldyn_conserv
17  !$OMP THREADPRIVATE(caldyn_conserv)
18
19  TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu
20
21CONTAINS
22
23  SUBROUTINE compute_planetvel(planetvel)
24    USE wind_mod
25    REAL(rstd),INTENT(OUT)  :: planetvel(iim*3*jjm)
26    REAL(rstd) :: ulon(iim*3*jjm)
27    REAL(rstd) :: ulat(iim*3*jjm)
28    REAL(rstd) :: lon,lat
29    INTEGER :: ij
30    DO ij=ij_begin_ext,ij_end_ext
31       ulon(ij+u_right)=radius*omega*cos(lat_e(ij+u_right))
32       ulat(ij+u_right)=0
33
34       ulon(ij+u_lup)=radius*omega*cos(lat_e(ij+u_lup))
35       ulat(ij+u_lup)=0
36
37       ulon(ij+u_ldown)=radius*omega*cos(lat_e(ij+u_ldown))
38       ulat(ij+u_ldown)=0
39    END DO
40    CALL compute_wind2D_perp_from_lonlat_compound(ulon, ulat, planetvel)
41  END SUBROUTINE compute_planetvel
42
43  !********************** compute_pvort = compute_theta + compute_pvort_only ******************
44
45  SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
46    USE icosa
47    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
48    USE exner_mod
49    USE trace
50    USE omp_para
51    IMPLICIT NONE
52    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
53    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
54    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
55    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
56    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
57    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
58    REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
59
60    INTEGER :: i,j,ij,l
61    REAL(rstd) :: etav,hv, m
62    !  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity 
63
64    CALL trace_start("compute_pvort") 
65
66    IF(caldyn_eta==eta_mass) THEN
67       CALL wait_message(req_ps)  ! COM00
68    ELSE
69       CALL wait_message(req_mass) ! COM00
70    END IF
71    CALL wait_message(req_theta_rhodz) ! COM01
72
73    IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
74       DO l = ll_begin,ll_end
75          CALL test_message(req_u) 
76          !DIR$ SIMD
77          DO ij=ij_begin_ext,ij_end_ext
78             m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
79             rhodz(ij,l) = m
80             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
81          ENDDO
82       ENDDO
83    ELSE ! Compute only theta
84       DO l = ll_begin,ll_end
85          CALL test_message(req_u) 
86          !DIR$ SIMD
87          DO ij=ij_begin_ext,ij_end_ext
88             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
89          ENDDO
90       ENDDO
91    END IF
92
93    CALL wait_message(req_u) ! COM02
94
95!!! Compute shallow-water potential vorticity
96    DO l = ll_begin,ll_end
97       !DIR$ SIMD
98       DO ij=ij_begin_ext,ij_end_ext
99          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
100               + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
101               - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
102
103          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
104               + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
105               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
106
107          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
108
109          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
110               + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
111               - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
112
113          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
114               + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
115               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
116
117          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
118
119       ENDDO
120
121       !DIR$ SIMD
122       DO ij=ij_begin,ij_end
123          qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l)) 
124          qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l)) 
125          qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l)) 
126       END DO
127
128    ENDDO
129
130    CALL trace_end("compute_pvort")
131  END SUBROUTINE compute_pvort
132
133  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
134    USE icosa
135    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
136    USE exner_mod
137    USE trace
138    USE omp_para
139    IMPLICIT NONE
140    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
141    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
142    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
143    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
144    INTEGER :: ij,l
145    REAL(rstd) :: m
146    CALL trace_start("compute_theta") 
147
148    IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
149       DO l = ll_begin,ll_end
150          !DIR$ SIMD
151          DO ij=ij_begin_ext,ij_end_ext
152             m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
153             rhodz(ij,l) = m
154             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
155          ENDDO
156       ENDDO
157    ELSE ! Compute only theta
158       DO l = ll_begin,ll_end
159          !DIR$ SIMD
160          DO ij=ij_begin_ext,ij_end_ext
161             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
162          ENDDO
163       ENDDO
164    END IF
165
166    CALL trace_end("compute_theta")
167  END SUBROUTINE compute_theta
168
169  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
170    USE icosa
171    USE exner_mod
172    USE trace
173    USE omp_para
174    IMPLICIT NONE
175    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
176    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
177    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
178    REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
179
180    INTEGER :: ij,l
181    REAL(rstd) :: etav,hv
182
183    CALL trace_start("compute_pvort_only") 
184!!! Compute shallow-water potential vorticity
185    DO l = ll_begin,ll_end
186       !DIR$ SIMD
187       DO ij=ij_begin_ext,ij_end_ext
188          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
189               + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
190               - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
191
192          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
193               + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
194               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
195
196          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
197
198          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
199               + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
200               - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
201
202          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
203               + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
204               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
205
206          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
207
208       ENDDO
209
210       !DIR$ SIMD
211       DO ij=ij_begin,ij_end
212          qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l)) 
213          qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l)) 
214          qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l)) 
215       END DO
216
217    ENDDO
218
219    CALL trace_end("compute_pvort_only")
220
221  END SUBROUTINE compute_pvort_only
222
223  !**************************** Geopotential *****************************
224
225  SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot) 
226    USE icosa
227    USE disvert_mod
228    USE exner_mod
229    USE trace
230    USE omp_para
231    IMPLICIT NONE
232    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
233    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
234    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature
235    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)       ! Exner function
236    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
237
238    INTEGER :: i,j,ij,l
239    REAL(rstd) :: p_ik, exner_ik
240    INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext
241
242
243    CALL trace_start("compute_geopot")
244
245    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
246    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1
247    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
248
249    IF(caldyn_eta==eta_mass) THEN
250
251!!! Compute exner function and geopotential
252       DO l = 1,llm
253          !DIR$ SIMD
254          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
255             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
256             !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
257             exner_ik = cpp * (p_ik/preff) ** kappa
258             pk(ij,l) = exner_ik
259             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
260             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
261          ENDDO
262       ENDDO
263       !       ENDIF
264    ELSE 
265       ! We are using a Lagrangian vertical coordinate
266       ! Pressure must be computed first top-down (temporarily stored in pk)
267       ! Then Exner pressure and geopotential are computed bottom-up
268       ! Notice that the computation below should work also when caldyn_eta=eta_mass         
269
270       IF(boussinesq) THEN ! compute only geopotential : pressure pk will be computed in compute_caldyn_horiz
271          ! specific volume 1 = dphi/g/rhodz
272          !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
273          DO l = 1,llm
274             !DIR$ SIMD
275             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
276                geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
277             ENDDO
278          ENDDO
279       ELSE ! non-Boussinesq, compute geopotential and Exner pressure
280          ! uppermost layer
281
282          !DIR$ SIMD
283          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
284             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
285          END DO
286          ! other layers
287          DO l = llm-1, 1, -1
288             !DIR$ SIMD
289             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
290                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
291             END DO
292          END DO
293          ! surface pressure (for diagnostics)
294          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
295             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
296          END DO
297
298          ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
299          DO l = 1,llm
300             !DIR$ SIMD
301             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
302                p_ik = pk(ij,l)
303                exner_ik = cpp * (p_ik/preff) ** kappa
304                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
305                pk(ij,l) = exner_ik
306             ENDDO
307          ENDDO
308       END IF
309
310    END IF
311
312    !ym flush geopot
313    !$OMP BARRIER
314
315    CALL trace_end("compute_geopot")
316
317  END SUBROUTINE compute_geopot
318
319  !************************* caldyn_horiz = caldyn_fast + caldyn_slow **********************
320
321  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du)
322    USE icosa
323    USE disvert_mod
324    USE exner_mod
325    USE trace
326    USE omp_para
327    IMPLICIT NONE
328    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
329    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
330    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
331    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
332    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
333    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
334
335    REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s
336    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
337    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
338    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
339
340    REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi)
341    REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
342    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
343    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
344
345    INTEGER :: i,j,ij,l
346    REAL(rstd) :: ww,uu
347
348    CALL trace_start("compute_caldyn_horiz")
349
350    !    CALL wait_message(req_theta_rhodz)
351
352    DO l = ll_begin, ll_end
353!!!  Compute mass and theta fluxes
354       IF (caldyn_conserv==energy) CALL test_message(req_qu) 
355       !DIR$ SIMD
356       DO ij=ij_begin_ext,ij_end_ext         
357          hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
358          hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
359          hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
360
361          Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
362          Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
363          Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
364       ENDDO
365
366!!! compute horizontal divergence of fluxes
367       !DIR$ SIMD
368       DO ij=ij_begin,ij_end         
369          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
370          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
371               ne_rup*hflux(ij+u_rup,l)       +  & 
372               ne_lup*hflux(ij+u_lup,l)       +  & 
373               ne_left*hflux(ij+u_left,l)     +  & 
374               ne_ldown*hflux(ij+u_ldown,l)   +  & 
375               ne_rdown*hflux(ij+u_rdown,l))   
376
377          ! signe ? attention d (rho theta dz)
378          ! dtheta_rhodz =  -div(flux.theta)
379          dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
380               ne_rup*Ftheta(ij+u_rup,l)        +  & 
381               ne_lup*Ftheta(ij+u_lup,l)        +  & 
382               ne_left*Ftheta(ij+u_left,l)      +  & 
383               ne_ldown*Ftheta(ij+u_ldown,l)    +  & 
384               ne_rdown*Ftheta(ij+u_rdown,l))
385       ENDDO
386
387    END DO
388
389!!! Compute potential vorticity (Coriolis) contribution to du
390
391    SELECT CASE(caldyn_conserv)
392    CASE(energy) ! energy-conserving TRiSK
393
394       CALL wait_message(req_qu) ! COM03
395
396       DO l=ll_begin,ll_end
397          !DIR$ SIMD
398          DO ij=ij_begin,ij_end         
399
400             uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
401                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
402                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
403                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
404                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
405                  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))+         &
406                  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))+         &
407                  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))+         &
408                  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))+             &
409                  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))
410             du(ij+u_right,l) = .5*uu/de(ij+u_right)
411
412             uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
413                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
414                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
415                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
416                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
417                  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)) +          &
418                  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)) +              &
419                  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)) +              &
420                  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)) +            &
421                  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))
422             du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
423
424
425             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
426                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
427                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
428                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
429                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
430                  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)) +          &
431                  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)) +        &
432                  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)) +      &
433                  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)) +      &
434                  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))
435             du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
436
437          ENDDO
438       ENDDO
439
440    CASE(enstrophy) ! enstrophy-conserving TRiSK
441
442       DO l=ll_begin,ll_end
443          !DIR$ SIMD
444          DO ij=ij_begin,ij_end         
445
446             uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
447                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
448                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
449                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
450                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
451                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
452                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
453                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
454                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
455                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
456             du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
457
458
459             uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
460                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
461                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
462                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
463                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
464                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
465                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
466                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
467                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
468                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
469             du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
470
471             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
472                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
473                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
474                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
475                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
476                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
477                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
478                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
479                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
480                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
481             du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
482
483          ENDDO
484       ENDDO
485
486    CASE DEFAULT
487       STOP
488    END SELECT
489
490!!! Compute bernouilli term = Kinetic Energy + geopotential
491    IF(boussinesq) THEN
492       ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)
493       ! uppermost layer
494       !DIR$ SIMD
495       DO ij=ij_begin_ext,ij_end_ext         
496          pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
497       END DO
498       ! other layers
499       DO l = llm-1, 1, -1
500          !          !$OMP DO SCHEDULE(STATIC)
501          !DIR$ SIMD
502          DO ij=ij_begin_ext,ij_end_ext         
503             pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
504          END DO
505       END DO
506       ! surface pressure (for diagnostics) FIXME
507       ! DO ij=ij_begin_ext,ij_end_ext         
508       !   ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1)
509       ! END DO
510       ! now pk contains the Lagrange multiplier (pressure)
511
512       DO l=ll_begin,ll_end
513          !DIR$ SIMD
514          DO ij=ij_begin,ij_end         
515
516             berni(ij,l) = pk(ij,l) + &
517                  1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
518                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
519                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
520                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
521                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
522                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
523             ! from now on pk contains the vertically-averaged geopotential
524             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
525          ENDDO
526       ENDDO
527
528    ELSE ! compressible
529
530       DO l=ll_begin,ll_end
531          !DIR$ SIMD
532          DO ij=ij_begin,ij_end         
533
534             berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
535                  + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
536                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
537                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
538                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
539                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
540                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
541          ENDDO
542       ENDDO
543
544    END IF ! Boussinesq/compressible
545
546!!! Add gradients of Bernoulli and Exner functions to du
547    DO l=ll_begin,ll_end
548       !DIR$ SIMD
549       DO ij=ij_begin,ij_end         
550
551          du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
552               0.5*(theta(ij,l)+theta(ij+t_right,l))                &
553               *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
554               + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
555
556
557          du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
558               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
559               *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
560               + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
561
562          du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
563               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
564               *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
565               + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
566
567       ENDDO
568    ENDDO
569
570    CALL trace_end("compute_caldyn_horiz")
571
572  END SUBROUTINE compute_caldyn_horiz
573
574  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot, du)
575    USE icosa
576    USE disvert_mod
577    USE exner_mod
578    USE trace
579    USE omp_para
580    IMPLICIT NONE
581    REAL(rstd), INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs
582    REAL(rstd),INTENT(INOUT)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
583    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
584    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
585    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
586    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
587    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
588    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
589
590    INTEGER :: i,j,ij,l
591
592    CALL trace_start("compute_caldyn_fast")
593
594    !    CALL wait_message(req_theta_rhodz)
595   
596    ! Compute bernouilli term
597    IF(boussinesq) THEN
598       ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)
599       ! uppermost layer
600       !DIR$ SIMD
601       DO ij=ij_begin_ext,ij_end_ext         
602          pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
603       END DO
604       ! other layers
605       DO l = llm-1, 1, -1
606          !          !$OMP DO SCHEDULE(STATIC)
607          !DIR$ SIMD
608          DO ij=ij_begin_ext,ij_end_ext         
609             pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
610          END DO
611       END DO
612       ! now pk contains the Lagrange multiplier (pressure)
613       DO l=ll_begin,ll_end
614          !DIR$ SIMD
615          DO ij=ij_begin,ij_end         
616             berni(ij,l) = pk(ij,l)
617             ! from now on pk contains the vertically-averaged geopotential
618             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
619          ENDDO
620       ENDDO
621
622    ELSE ! compressible
623
624       DO l=ll_begin,ll_end
625          !DIR$ SIMD
626          DO ij=ij_begin,ij_end         
627             berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
628          ENDDO
629       ENDDO
630
631    END IF ! Boussinesq/compressible
632
633!!! u:=u+tau*du, du = gradients of Bernoulli and Exner functions
634    DO l=ll_begin,ll_end
635       !DIR$ SIMD
636       DO ij=ij_begin,ij_end         
637
638          du(ij+u_right,l) = 1/de(ij+u_right) * (             &
639               0.5*(theta(ij,l)+theta(ij+t_right,l))                &
640               *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
641               + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
642          u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
643
644
645          du(ij+u_lup,l) = 1/de(ij+u_lup) * (                  &
646               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
647               *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
648               + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
649          u(ij+u_lup,l) = u(ij+u_lup,l) + tau*du(ij+u_lup,l)
650
651          du(ij+u_ldown,l) = 1/de(ij+u_ldown) * (             &
652               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
653               *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
654               + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
655          u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
656
657       ENDDO
658    ENDDO
659
660    CALL trace_end("compute_caldyn_fast")
661
662  END SUBROUTINE compute_caldyn_fast
663
664  SUBROUTINE compute_caldyn_slow(u,rhodz,qu,theta, hflux,convm, dtheta_rhodz, du)
665    USE icosa
666    USE disvert_mod
667    USE exner_mod
668    USE trace
669    USE omp_para
670    IMPLICIT NONE
671    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
672    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
673    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
674    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
675
676    REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s
677    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
678    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
679    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
680
681    REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi)
682    REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
683    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
684    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
685
686    INTEGER :: i,j,ij,l
687    REAL(rstd) :: ww,uu
688
689    CALL trace_start("compute_caldyn_slow")
690
691    !    CALL wait_message(req_theta_rhodz)
692
693    DO l = ll_begin, ll_end
694!!!  Compute mass and theta fluxes
695       IF (caldyn_conserv==energy) CALL test_message(req_qu) 
696       !DIR$ SIMD
697       DO ij=ij_begin_ext,ij_end_ext         
698          hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
699          hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
700          hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
701
702          Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
703          Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
704          Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
705       ENDDO
706
707!!! compute horizontal divergence of fluxes
708       !DIR$ SIMD
709       DO ij=ij_begin,ij_end         
710          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
711          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
712               ne_rup*hflux(ij+u_rup,l)       +  & 
713               ne_lup*hflux(ij+u_lup,l)       +  & 
714               ne_left*hflux(ij+u_left,l)     +  & 
715               ne_ldown*hflux(ij+u_ldown,l)   +  & 
716               ne_rdown*hflux(ij+u_rdown,l))   
717
718          ! signe ? attention d (rho theta dz)
719          ! dtheta_rhodz =  -div(flux.theta)
720          dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
721               ne_rup*Ftheta(ij+u_rup,l)        +  & 
722               ne_lup*Ftheta(ij+u_lup,l)        +  & 
723               ne_left*Ftheta(ij+u_left,l)      +  & 
724               ne_ldown*Ftheta(ij+u_ldown,l)    +  & 
725               ne_rdown*Ftheta(ij+u_rdown,l))
726       ENDDO
727
728    END DO
729
730!!! Compute potential vorticity (Coriolis) contribution to du
731
732    SELECT CASE(caldyn_conserv)
733    CASE(energy) ! energy-conserving TRiSK
734
735       CALL wait_message(req_qu)
736
737       DO l=ll_begin,ll_end
738          !DIR$ SIMD
739          DO ij=ij_begin,ij_end         
740
741             uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
742                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
743                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
744                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
745                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
746                  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))+         &
747                  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))+         &
748                  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))+         &
749                  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))+             &
750                  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))
751             du(ij+u_right,l) = .5*uu/de(ij+u_right)
752
753             uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
754                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
755                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
756                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
757                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
758                  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)) +          &
759                  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)) +              &
760                  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)) +              &
761                  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)) +            &
762                  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))
763             du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
764
765
766             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
767                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
768                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
769                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
770                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
771                  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)) +          &
772                  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)) +        &
773                  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)) +      &
774                  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)) +      &
775                  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))
776             du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
777
778          ENDDO
779       ENDDO
780
781    CASE(enstrophy) ! enstrophy-conserving TRiSK
782
783       DO l=ll_begin,ll_end
784          !DIR$ SIMD
785          DO ij=ij_begin,ij_end         
786
787             uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
788                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
789                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
790                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
791                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
792                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
793                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
794                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
795                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
796                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
797             du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
798
799
800             uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
801                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
802                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
803                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
804                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
805                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
806                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
807                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
808                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
809                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
810             du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
811
812             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
813                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
814                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
815                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
816                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
817                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
818                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
819                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
820                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
821                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
822             du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
823
824          ENDDO
825       ENDDO
826
827    CASE DEFAULT
828       STOP
829    END SELECT
830
831    ! Compute bernouilli term = Kinetic Energy
832    IF(boussinesq) THEN
833
834       DO l=ll_begin,ll_end
835          !DIR$ SIMD
836          DO ij=ij_begin,ij_end         
837
838             berni(ij,l) = &
839                  1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
840                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
841                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
842                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
843                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
844                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
845          ENDDO
846       ENDDO
847
848    ELSE ! compressible
849
850       DO l=ll_begin,ll_end
851          !DIR$ SIMD
852          DO ij=ij_begin,ij_end         
853
854             berni(ij,l) =  &
855                  + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
856                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
857                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
858                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
859                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
860                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
861          ENDDO
862       ENDDO
863
864    END IF ! Boussinesq/compressible
865
866!!! Add gradients of Bernoulli and Exner functions to du
867    DO l=ll_begin,ll_end
868       !DIR$ SIMD
869       DO ij=ij_begin,ij_end
870          du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right)            &
871               * ( ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
872          du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup)                 &
873               * ( ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
874          du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown)            &
875               * ( ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
876       ENDDO
877    ENDDO
878
879    CALL trace_end("compute_caldyn_slow")
880
881  END SUBROUTINE compute_caldyn_slow
882
883  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
884    USE icosa
885    USE disvert_mod
886    USE exner_mod
887    USE trace
888    USE omp_para
889    IMPLICIT NONE
890    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
891    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
892    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
893    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
894
895    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
896    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
897    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
898    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm)
899    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
900
901    ! temporary variable   
902    INTEGER :: i,j,ij,l
903    REAL(rstd) :: p_ik, exner_ik
904    INTEGER    :: ij_omp_begin, ij_omp_end
905
906
907    CALL trace_start("compute_caldyn_vert")
908
909    CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin,ij_omp_end)
910    ij_omp_begin=ij_omp_begin+ij_begin-1
911    ij_omp_end=ij_omp_end+ij_begin-1
912
913    !    REAL(rstd) :: wwuu(iim*3*jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp
914    ! need to be understood
915
916    !    wwuu=wwuu_out
917    CALL trace_start("compute_caldyn_vert")
918
919    !$OMP BARRIER   
920!!! cumulate mass flux convergence from top to bottom
921    !  IF (is_omp_level_master) THEN
922    DO  l = llm-1, 1, -1
923       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)
924
925!!$OMP DO SCHEDULE(STATIC)
926       !DIR$ SIMD
927       DO ij=ij_omp_begin,ij_omp_end         
928          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
929       ENDDO
930    ENDDO
931    !  ENDIF
932
933    !$OMP BARRIER
934    ! FLUSH on convm
935!!!!!!!!!!!!!!!!!!!!!!!!!
936
937    ! compute dps
938    IF (is_omp_first_level) THEN
939       !DIR$ SIMD
940       DO ij=ij_begin,ij_end         
941          ! dps/dt = -int(div flux)dz
942          dps(ij) = convm(ij,1) * g 
943       ENDDO
944    ENDIF
945
946!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
947    DO l=ll_beginp1,ll_end
948       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)
949       !DIR$ SIMD
950       DO ij=ij_begin,ij_end         
951          ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
952          ! => w>0 for upward transport
953          wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
954       ENDDO
955    ENDDO
956
957    !--> flush wflux 
958    !$OMP BARRIER
959
960    DO l=ll_begin,ll_endm1
961       !DIR$ SIMD
962       DO ij=ij_begin,ij_end         
963          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
964       ENDDO
965    ENDDO
966
967    DO l=ll_beginp1,ll_end
968       !DIR$ SIMD
969       DO ij=ij_begin,ij_end         
970          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
971       ENDDO
972    ENDDO
973
974
975    ! Compute vertical transport
976    DO l=ll_beginp1,ll_end
977       !DIR$ SIMD
978       DO ij=ij_begin,ij_end         
979          wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1))
980          wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
981          wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
982       ENDDO
983    ENDDO
984
985    !--> flush wwuu 
986    !$OMP BARRIER
987
988    ! Add vertical transport to du
989    DO l=ll_begin,ll_end
990       !DIR$ SIMD
991       DO ij=ij_begin,ij_end         
992          du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
993          du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
994          du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
995       ENDDO
996    ENDDO
997
998    !  DO l=ll_beginp1,ll_end
999    !!DIR$ SIMD
1000    !      DO ij=ij_begin,ij_end         
1001    !        wwuu_out(ij+u_right,l) = wwuu(ij+u_right,l)
1002    !        wwuu_out(ij+u_lup,l)   = wwuu(ij+u_lup,l)
1003    !        wwuu_out(ij+u_ldown,l) = wwuu(ij+u_ldown,l)
1004    !     ENDDO
1005    !   ENDDO
1006
1007    CALL trace_end("compute_caldyn_vert")
1008
1009  END SUBROUTINE compute_caldyn_vert
1010
1011END MODULE caldyn_kernels_mod
Note: See TracBrowser for help on using the repository browser.