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

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

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