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

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

devel : separate modules for caldyn_vert and caldyn_vert_NH
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 : id_coriolis, primal_num, dual_num, edge_num, &
66          left, right,primal_deg,primal_edge,primal_ne,trisk_deg,wee,trisk,Ai, &
67          enter_trace, exit_trace
68    FIELD_U     :: hflux, Ftheta, qu, du ! IN, BUF, IN, INOUT
69    FIELD_MASS  :: convm                 ! BUF
70    FIELD_THETA :: theta, dtheta_rhodz   ! IN, OUT
71    DECLARE_INDICES
72    DECLARE_EDGES
73    NUM :: divF, du_trisk
74    START_TRACE(id_coriolis, 3,4,0) ! primal, dual, edge
75#include "../kernels_unst/coriolis.k90"
76    STOP_TRACE
77  END SUBROUTINE compute_caldyn_coriolis_unst
78
79  SUBROUTINE compute_caldyn_Coriolis_hex(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
80    USE icosa
81    USE caldyn_vars_mod
82    REAL(rstd),INTENT(IN)    :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
83    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars
84    REAL(rstd),INTENT(IN)    :: qu(3*iim*jjm,llm)
85    REAL(rstd), INTENT(OUT)  :: Ftheta(3*iim*jjm,llm)  ! potential temperature flux
86    REAL(rstd),INTENT(OUT)   :: convm(iim*jjm,llm)  ! mass flux convergence
87    REAL(rstd),INTENT(OUT)   :: dtheta_rhodz(iim*jjm,llm,nqdyn)
88    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
89    REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
90    INTEGER :: ij,iq,l,kdown
91    CALL trace_start("compute_caldyn_Coriolis")
92
93#include "../kernels_hex/coriolis.k90"
94
95    CALL trace_end("compute_caldyn_Coriolis")
96  END SUBROUTINE compute_caldyn_Coriolis_hex
97
98  SUBROUTINE compute_caldyn_Coriolis_manual(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
99    USE icosa
100    USE caldyn_vars_mod
101    REAL(rstd),INTENT(IN)    :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
102    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars
103    REAL(rstd),INTENT(IN)    :: qu(3*iim*jjm,llm)
104    REAL(rstd),INTENT(IN)    :: Ftheta(3*iim*jjm,llm)  ! ignored in favor of local buffer
105    REAL(rstd),INTENT(OUT)   :: convm(iim*jjm,llm)  ! mass flux convergence
106    REAL(rstd),INTENT(OUT)   :: dtheta_rhodz(iim*jjm,llm,nqdyn)
107    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
108    REAL(rstd)  :: buf_Ftheta(3*iim*jjm)  ! local buffer for potential temperature flux
109    REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
110    INTEGER :: ij,iq,l,kdown
111
112    CALL trace_start("compute_caldyn_Coriolis")
113
114#define FTHETA(ij) buf_Ftheta(ij)
115
116    DO l=ll_begin, ll_end
117       ! compute theta flux
118       DO iq=1,nqdyn
119       !DIR$ SIMD
120          DO ij=ij_begin_ext,ij_end_ext     
121             FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
122                                  * hflux(ij+u_right,l)
123             FTHETA(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
124                  * hflux(ij+u_lup,l)
125             FTHETA(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
126                  * hflux(ij+u_ldown,l)
127          END DO
128          ! horizontal divergence of fluxes
129       !DIR$ SIMD
130          DO ij=ij_begin,ij_end         
131             ! dtheta_rhodz =  -div(flux.theta)
132             dtheta_rhodz(ij,l,iq)= &
133                  -1./Ai(ij)*(ne_right*FTHETA(ij+u_right)    +  &
134                  ne_rup*FTHETA(ij+u_rup)        +  & 
135                  ne_lup*FTHETA(ij+u_lup)        +  & 
136                  ne_left*FTHETA(ij+u_left)      +  & 
137                  ne_ldown*FTHETA(ij+u_ldown)    +  &
138                  ne_rdown*FTHETA(ij+u_rdown) )
139          END DO
140       END DO
141
142       !DIR$ SIMD
143       DO ij=ij_begin,ij_end         
144          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
145          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
146               ne_rup*hflux(ij+u_rup,l)       +  & 
147               ne_lup*hflux(ij+u_lup,l)       +  & 
148               ne_left*hflux(ij+u_left,l)     +  & 
149               ne_ldown*hflux(ij+u_ldown,l)   +  & 
150               ne_rdown*hflux(ij+u_rdown,l))
151       END DO ! ij
152    END DO ! llm
153
154!!! Compute potential vorticity (Coriolis) contribution to du
155    SELECT CASE(caldyn_conserv)
156
157    CASE(conserv_energy) ! energy-conserving TRiSK
158
159       DO l=ll_begin,ll_end
160          !DIR$ SIMD
161          DO ij=ij_begin,ij_end         
162             uu_right = &
163                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
164                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
165                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
166                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
167                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
168                  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))+         &
169                  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))+         &
170                  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))+         &
171                  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))+             &
172                  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))
173             uu_lup = &
174                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
175                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
176                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
177                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
178                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
179                  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)) +          &
180                  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)) +              &
181                  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)) +              &
182                  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)) +            &
183                  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))
184             uu_ldown = &
185                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
186                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
187                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
188                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
189                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
190                  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)) +          &
191                  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)) +        &
192                  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)) +      &
193                  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)) +      &
194                  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))
195             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right
196             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup
197             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown
198          ENDDO
199       ENDDO
200
201    CASE(conserv_gassmann) ! energy-conserving TRiSK modified by Gassmann (2018)
202
203       DO l=ll_begin,ll_end
204          !DIR$ SIMD
205          DO ij=ij_begin,ij_end         
206             uu_right = &
207                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)  *qu(ij+t_right+u_lup,l)+                 &
208                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)  *qu(ij+u_rup,l)+                         &
209                .5*wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+     &
210                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*qu(ij+u_rdown,l)+                       &
211                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*qu(ij+t_right+u_ldown,l)+               &
212                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*qu(ij+u_rdown,l)+               &
213                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*qu(ij+t_right+u_ldown,l)+       &
214               .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))+         &
215                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*qu(ij+t_right+u_lup,l)+           &
216                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*qu(ij+u_rup,l)
217             uu_lup = &
218                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*qu(ij+t_lup+u_ldown,l) +                   &
219                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*qu(ij+u_left,l) +                         &
220               .5*wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +       &
221                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*qu(ij+u_rup,l) +                          &
222                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*qu(ij+t_lup+u_right,l)+                     & 
223                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*qu(ij+u_rup,l) +                   &
224                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*qu(ij+t_lup+u_right,l) +              &
225               .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)) + &
226                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*qu(ij+t_lup+u_ldown,l) +            &
227                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*qu(ij+u_left,l)
228             uu_ldown = &
229                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*qu(ij+t_ldown,l+u_right) +               &
230                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*qu(ij+u_rdown,l) +                       &
231               .5*wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +        &
232                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*qu(ij+u_left,l) +                          &
233                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*qu(ij+t_ldown+u_lup,l) +                  & 
234                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*qu(ij+u_left,l) +                    &
235                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*qu(ij+t_ldown+u_lup,l) +          &
236               .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)) +      &
237                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*qu(ij+t_ldown+u_right,l) +      &
238                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*qu(ij+u_rdown,l)
239             du(ij+u_right,l) = du(ij+u_right,l) + uu_right
240             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup
241             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + uu_ldown
242          ENDDO
243       ENDDO
244
245    CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK
246
247       DO l=ll_begin,ll_end
248          !DIR$ SIMD
249          DO ij=ij_begin,ij_end         
250             uu_right = &
251                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
252                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
253                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
254                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
255                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
256                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
257                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
258                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
259                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
260                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
261             uu_lup = &
262                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
263                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
264                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
265                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
266                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
267                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
268                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
269                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
270                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
271                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
272             uu_ldown = &
273                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
274                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
275                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
276                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
277                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
278                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
279                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
280                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
281                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
282                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
283
284             du(ij+u_right,l) = du(ij+u_right,l) + uu_right*qu(ij+u_right,l)
285             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup*qu(ij+u_lup,l)     
286             du(ij+u_ldown,l) = du(ij+u_ldown,l) + uu_ldown*qu(ij+u_ldown,l) 
287          END DO
288       END DO
289    CASE DEFAULT
290       STOP
291    END SELECT
292#undef FTHETA
293
294    CALL trace_end("compute_caldyn_Coriolis")
295
296  END SUBROUTINE compute_caldyn_Coriolis_manual
297 
298END MODULE compute_caldyn_Coriolis_mod
Note: See TracBrowser for help on using the repository browser.