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

Last change on this file since 876 was 876, checked in by jisesh, 5 years ago

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