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

Last change on this file since 428 was 428, checked in by dubos, 8 years ago

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