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

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

devel: separate module for compute_caldyn_Coriolis
File size: 12.9 KB
Line 
1MODULE compute_caldyn_Coriolis_mod
2  USE grid_param, ONLY : llm
3  IMPLICIT NONE
4  PRIVATE
5
6  PUBLIC :: compute_caldyn_Coriolis
7
8CONTAINS
9
10  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
11    USE icosa
12    USE trace
13    USE caldyn_vars_mod
14    USE omp_para, ONLY : ll_begin, ll_end
15    REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
16    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) ! active scalars
17    REAL(rstd),INTENT(IN)  :: qu(3*iim*jjm,llm)
18    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
19    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
20    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
21   
22    REAL(rstd) :: Ftheta(3*iim*jjm,llm)  ! potential temperature flux
23    REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
24    INTEGER :: ij,iq,l,kdown
25
26    CALL trace_start("compute_caldyn_Coriolis")
27
28    IF(dysl_caldyn_coriolis) THEN
29
30#include "../kernels_hex/coriolis.k90"
31
32    ELSE
33#define FTHETA(ij) Ftheta(ij,1)
34
35    DO l=ll_begin, ll_end
36       ! compute theta flux
37       DO iq=1,nqdyn
38       !DIR$ SIMD
39          DO ij=ij_begin_ext,ij_end_ext     
40             FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
41                                  * hflux(ij+u_right,l)
42             FTHETA(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
43                  * hflux(ij+u_lup,l)
44             FTHETA(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
45                  * hflux(ij+u_ldown,l)
46          END DO
47          ! horizontal divergence of fluxes
48       !DIR$ SIMD
49          DO ij=ij_begin,ij_end         
50             ! dtheta_rhodz =  -div(flux.theta)
51             dtheta_rhodz(ij,l,iq)= &
52                  -1./Ai(ij)*(ne_right*FTHETA(ij+u_right)    +  &
53                  ne_rup*FTHETA(ij+u_rup)        +  & 
54                  ne_lup*FTHETA(ij+u_lup)        +  & 
55                  ne_left*FTHETA(ij+u_left)      +  & 
56                  ne_ldown*FTHETA(ij+u_ldown)    +  &
57                  ne_rdown*FTHETA(ij+u_rdown) )
58          END DO
59       END DO
60
61       !DIR$ SIMD
62       DO ij=ij_begin,ij_end         
63          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
64          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
65               ne_rup*hflux(ij+u_rup,l)       +  & 
66               ne_lup*hflux(ij+u_lup,l)       +  & 
67               ne_left*hflux(ij+u_left,l)     +  & 
68               ne_ldown*hflux(ij+u_ldown,l)   +  & 
69               ne_rdown*hflux(ij+u_rdown,l))
70       END DO ! ij
71    END DO ! llm
72
73!!! Compute potential vorticity (Coriolis) contribution to du
74    SELECT CASE(caldyn_conserv)
75
76    CASE(conserv_energy) ! energy-conserving TRiSK
77
78       DO l=ll_begin,ll_end
79          !DIR$ SIMD
80          DO ij=ij_begin,ij_end         
81             uu_right = &
82                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
83                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
84                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
85                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
86                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
87                  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))+         &
88                  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))+         &
89                  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))+         &
90                  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))+             &
91                  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))
92             uu_lup = &
93                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
94                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
95                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
96                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
97                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
98                  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)) +          &
99                  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)) +              &
100                  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)) +              &
101                  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)) +            &
102                  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))
103             uu_ldown = &
104                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
105                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
106                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
107                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
108                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
109                  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)) +          &
110                  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)) +        &
111                  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)) +      &
112                  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)) +      &
113                  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))
114             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right
115             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup
116             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown
117          ENDDO
118       ENDDO
119
120    CASE(conserv_gassmann) ! energy-conserving TRiSK modified by Gassmann (2018)
121
122       DO l=ll_begin,ll_end
123          !DIR$ SIMD
124          DO ij=ij_begin,ij_end         
125             uu_right = &
126                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)  *qu(ij+t_right+u_lup,l)+                 &
127                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)  *qu(ij+u_rup,l)+                         &
128                .5*wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+     &
129                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*qu(ij+u_rdown,l)+                       &
130                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*qu(ij+t_right+u_ldown,l)+               &
131                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*qu(ij+u_rdown,l)+               &
132                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*qu(ij+t_right+u_ldown,l)+       &
133               .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))+         &
134                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*qu(ij+t_right+u_lup,l)+           &
135                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*qu(ij+u_rup,l)
136             uu_lup = &
137                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*qu(ij+t_lup+u_ldown,l) +                   &
138                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*qu(ij+u_left,l) +                         &
139               .5*wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +       &
140                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*qu(ij+u_rup,l) +                          &
141                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*qu(ij+t_lup+u_right,l)+                     & 
142                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*qu(ij+u_rup,l) +                   &
143                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*qu(ij+t_lup+u_right,l) +              &
144               .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)) + &
145                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*qu(ij+t_lup+u_ldown,l) +            &
146                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*qu(ij+u_left,l)
147             uu_ldown = &
148                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*qu(ij+t_ldown,l+u_right) +               &
149                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*qu(ij+u_rdown,l) +                       &
150               .5*wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +        &
151                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*qu(ij+u_left,l) +                          &
152                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*qu(ij+t_ldown+u_lup,l) +                  & 
153                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*qu(ij+u_left,l) +                    &
154                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*qu(ij+t_ldown+u_lup,l) +          &
155               .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)) +      &
156                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*qu(ij+t_ldown+u_right,l) +      &
157                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*qu(ij+u_rdown,l)
158             du(ij+u_right,l) = du(ij+u_right,l) + uu_right
159             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup
160             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + uu_ldown
161          ENDDO
162       ENDDO
163
164    CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK
165
166       DO l=ll_begin,ll_end
167          !DIR$ SIMD
168          DO ij=ij_begin,ij_end         
169             uu_right = &
170                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
171                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
172                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
173                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
174                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
175                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
176                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
177                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
178                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
179                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
180             uu_lup = &
181                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
182                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
183                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
184                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
185                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
186                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
187                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
188                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
189                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
190                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
191             uu_ldown = &
192                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
193                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
194                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
195                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
196                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
197                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
198                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
199                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
200                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
201                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
202
203             du(ij+u_right,l) = du(ij+u_right,l) + uu_right*qu(ij+u_right,l)
204             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup*qu(ij+u_lup,l)     
205             du(ij+u_ldown,l) = du(ij+u_ldown,l) + uu_ldown*qu(ij+u_ldown,l) 
206          END DO
207       END DO
208    CASE DEFAULT
209       STOP
210    END SELECT
211#undef FTHETA
212
213    END IF ! dysl
214
215    CALL trace_end("compute_caldyn_Coriolis")
216
217  END SUBROUTINE compute_caldyn_Coriolis
218 
219END MODULE compute_caldyn_Coriolis_mod
Note: See TracBrowser for help on using the repository browser.