source: codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90 @ 939

Last change on this file since 939 was 939, checked in by dubos, 5 years ago

devel : cleanup USE data_unstructured_mod
File size: 15.5 KB
Line 
1MODULE compute_caldyn_Coriolis_mod
2  USE prec, ONLY : rstd
3  USE grid_param
4  USE earth_const
5  USE disvert_mod
6  USE omp_para
7  USE trace
8  IMPLICIT NONE
9  PRIVATE
10
11#include "../unstructured/unstructured.h90"
12
13  PUBLIC :: compute_caldyn_Coriolis_manual, &
14       compute_caldyn_Coriolis_unst, compute_caldyn_Coriolis_hex
15
16CONTAINS
17
18#ifdef BEGIN_DYSL
19
20KERNEL(coriolis)
21!
22  DO iq=1,nqdyn
23    FORALL_CELLS_EXT()
24      ON_EDGES
25        Ftheta(EDGE) = .5*(theta(CELL1,iq)+theta(CELL2,iq))*hflux(EDGE)
26      END_BLOCK
27    END_BLOCK
28    FORALL_CELLS()
29      ON_PRIMAL
30        divF=0.
31        FORALL_EDGES
32          divF = divF + Ftheta(EDGE)*SIGN
33        END_BLOCK
34        dtheta_rhodz(CELL,iq) = -divF / AI
35      END_BLOCK
36    END_BLOCK
37  END DO ! iq
38!
39  FORALL_CELLS()
40    ON_PRIMAL
41      divF=0.
42      FORALL_EDGES
43        divF = divF + hflux(EDGE)*SIGN
44      END_BLOCK
45      convm(CELL) = -divF / AI
46    END_BLOCK
47  END_BLOCK
48!
49  FORALL_CELLS()
50    ON_EDGES
51      du_trisk=0.
52      FORALL_TRISK
53        du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK))
54      END_BLOCK
55      du(EDGE) = du(EDGE) + .5*du_trisk
56    END_BLOCK
57  END_BLOCK
58
59END_BLOCK
60
61#endif END_DYSL
62
63  SUBROUTINE compute_caldyn_coriolis_unst(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
64    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
65    USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
66         id_coriolis, left, right, primal_deg, primal_edge, primal_ne, &
67         trisk_deg, trisk, wee ! FIXME wee
68         
69    FIELD_U     :: hflux, Ftheta, qu, du ! IN, BUF, IN, INOUT
70    FIELD_MASS  :: convm                 ! BUF
71    FIELD_THETA :: theta, dtheta_rhodz   ! IN, OUT
72    DECLARE_INDICES
73    DECLARE_EDGES
74    NUM :: divF, du_trisk
75    START_TRACE(id_coriolis, 3,4,0) ! primal, dual, edge
76#include "../kernels_unst/coriolis.k90"
77    STOP_TRACE
78  END SUBROUTINE compute_caldyn_coriolis_unst
79
80  SUBROUTINE compute_caldyn_Coriolis_hex(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
81    USE icosa
82    USE caldyn_vars_mod
83    REAL(rstd),INTENT(IN)    :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
84    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars
85    REAL(rstd),INTENT(IN)    :: qu(3*iim*jjm,llm)
86    REAL(rstd), INTENT(OUT)  :: Ftheta(3*iim*jjm,llm)  ! potential temperature flux
87    REAL(rstd),INTENT(OUT)   :: convm(iim*jjm,llm)  ! mass flux convergence
88    REAL(rstd),INTENT(OUT)   :: dtheta_rhodz(iim*jjm,llm,nqdyn)
89    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
90    REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
91    INTEGER :: ij,iq,l,kdown
92    CALL trace_start("compute_caldyn_Coriolis")
93
94#include "../kernels_hex/coriolis.k90"
95
96    CALL trace_end("compute_caldyn_Coriolis")
97  END SUBROUTINE compute_caldyn_Coriolis_hex
98
99  SUBROUTINE compute_caldyn_Coriolis_manual(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
100    USE icosa
101    USE caldyn_vars_mod
102    REAL(rstd),INTENT(IN)    :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
103    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars
104    REAL(rstd),INTENT(IN)    :: qu(3*iim*jjm,llm)
105    REAL(rstd),INTENT(IN)    :: Ftheta(3*iim*jjm,llm)  ! ignored in favor of local buffer
106    REAL(rstd),INTENT(OUT)   :: convm(iim*jjm,llm)  ! mass flux convergence
107    REAL(rstd),INTENT(OUT)   :: dtheta_rhodz(iim*jjm,llm,nqdyn)
108    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
109    REAL(rstd)  :: buf_Ftheta(3*iim*jjm)  ! local buffer for potential temperature flux
110    REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
111    INTEGER :: ij,iq,l,kdown
112
113    CALL trace_start("compute_caldyn_Coriolis")
114
115#define FTHETA(ij) buf_Ftheta(ij)
116
117    DO l=ll_begin, ll_end
118       ! compute theta flux
119       DO iq=1,nqdyn
120       !DIR$ SIMD
121          DO ij=ij_begin_ext,ij_end_ext     
122             FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
123                                  * hflux(ij+u_right,l)
124             FTHETA(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
125                  * hflux(ij+u_lup,l)
126             FTHETA(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
127                  * hflux(ij+u_ldown,l)
128          END DO
129          ! horizontal divergence of fluxes
130       !DIR$ SIMD
131          DO ij=ij_begin,ij_end         
132             ! dtheta_rhodz =  -div(flux.theta)
133             dtheta_rhodz(ij,l,iq)= &
134                  -1./Ai(ij)*(ne_right*FTHETA(ij+u_right)    +  &
135                  ne_rup*FTHETA(ij+u_rup)        +  & 
136                  ne_lup*FTHETA(ij+u_lup)        +  & 
137                  ne_left*FTHETA(ij+u_left)      +  & 
138                  ne_ldown*FTHETA(ij+u_ldown)    +  &
139                  ne_rdown*FTHETA(ij+u_rdown) )
140          END DO
141       END DO
142
143       !DIR$ SIMD
144       DO ij=ij_begin,ij_end         
145          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
146          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
147               ne_rup*hflux(ij+u_rup,l)       +  & 
148               ne_lup*hflux(ij+u_lup,l)       +  & 
149               ne_left*hflux(ij+u_left,l)     +  & 
150               ne_ldown*hflux(ij+u_ldown,l)   +  & 
151               ne_rdown*hflux(ij+u_rdown,l))
152       END DO ! ij
153    END DO ! llm
154
155!!! Compute potential vorticity (Coriolis) contribution to du
156    SELECT CASE(caldyn_conserv)
157
158    CASE(conserv_energy) ! energy-conserving TRiSK
159
160       DO l=ll_begin,ll_end
161          !DIR$ SIMD
162          DO ij=ij_begin,ij_end         
163             uu_right = &
164                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
165                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
166                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
167                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
168                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
169                  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))+         &
170                  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))+         &
171                  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))+         &
172                  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))+             &
173                  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))
174             uu_lup = &
175                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
176                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
177                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
178                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
179                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
180                  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)) +          &
181                  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)) +              &
182                  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)) +              &
183                  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)) +            &
184                  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))
185             uu_ldown = &
186                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
187                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
188                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
189                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
190                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
191                  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)) +          &
192                  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)) +        &
193                  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)) +      &
194                  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)) +      &
195                  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))
196             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right
197             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup
198             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown
199          ENDDO
200       ENDDO
201
202    CASE(conserv_gassmann) ! energy-conserving TRiSK modified by Gassmann (2018)
203
204       DO l=ll_begin,ll_end
205          !DIR$ SIMD
206          DO ij=ij_begin,ij_end         
207             uu_right = &
208                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)  *qu(ij+t_right+u_lup,l)+                 &
209                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)  *qu(ij+u_rup,l)+                         &
210                .5*wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+     &
211                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*qu(ij+u_rdown,l)+                       &
212                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*qu(ij+t_right+u_ldown,l)+               &
213                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*qu(ij+u_rdown,l)+               &
214                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*qu(ij+t_right+u_ldown,l)+       &
215               .5*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))+         &
216                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*qu(ij+t_right+u_lup,l)+           &
217                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*qu(ij+u_rup,l)
218             uu_lup = &
219                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*qu(ij+t_lup+u_ldown,l) +                   &
220                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*qu(ij+u_left,l) +                         &
221               .5*wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +       &
222                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*qu(ij+u_rup,l) +                          &
223                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*qu(ij+t_lup+u_right,l)+                     & 
224                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*qu(ij+u_rup,l) +                   &
225                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*qu(ij+t_lup+u_right,l) +              &
226               .5*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)) + &
227                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*qu(ij+t_lup+u_ldown,l) +            &
228                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*qu(ij+u_left,l)
229             uu_ldown = &
230                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*qu(ij+t_ldown,l+u_right) +               &
231                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*qu(ij+u_rdown,l) +                       &
232               .5*wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +        &
233                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*qu(ij+u_left,l) +                          &
234                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*qu(ij+t_ldown+u_lup,l) +                  & 
235                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*qu(ij+u_left,l) +                    &
236                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*qu(ij+t_ldown+u_lup,l) +          &
237               .5*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)) +      &
238                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*qu(ij+t_ldown+u_right,l) +      &
239                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*qu(ij+u_rdown,l)
240             du(ij+u_right,l) = du(ij+u_right,l) + uu_right
241             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup
242             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + uu_ldown
243          ENDDO
244       ENDDO
245
246    CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK
247
248       DO l=ll_begin,ll_end
249          !DIR$ SIMD
250          DO ij=ij_begin,ij_end         
251             uu_right = &
252                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
253                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
254                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
255                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
256                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
257                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
258                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
259                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
260                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
261                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
262             uu_lup = &
263                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
264                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
265                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
266                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
267                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
268                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
269                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
270                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
271                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
272                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
273             uu_ldown = &
274                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
275                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
276                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
277                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
278                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
279                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
280                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
281                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
282                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
283                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
284
285             du(ij+u_right,l) = du(ij+u_right,l) + uu_right*qu(ij+u_right,l)
286             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup*qu(ij+u_lup,l)     
287             du(ij+u_ldown,l) = du(ij+u_ldown,l) + uu_ldown*qu(ij+u_ldown,l) 
288          END DO
289       END DO
290    CASE DEFAULT
291       STOP
292    END SELECT
293#undef FTHETA
294
295    CALL trace_end("compute_caldyn_Coriolis")
296
297  END SUBROUTINE compute_caldyn_Coriolis_manual
298 
299END MODULE compute_caldyn_Coriolis_mod
Note: See TracBrowser for help on using the repository browser.