source: codes/icosagcm/trunk/src/dynamics/caldyn_kernels.f90

Last change on this file was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

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