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

Last change on this file since 519 was 519, checked in by dubos, 7 years ago

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