source: codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90 @ 921

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

devel : unique interface for compute_caldyn_fast_X, compute_caldyn_slow_hydro_X
File size: 8.5 KB
Line 
1MODULE compute_caldyn_fast_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_fast_unst, compute_caldyn_fast_hex, compute_caldyn_fast_manual
14
15CONTAINS
16
17#ifdef BEGIN_DYSL
18
19KERNEL(caldyn_fast)
20!
21  SELECT CASE(caldyn_thermo)
22  CASE(thermo_boussinesq)
23    FORALL_CELLS()
24      ON_PRIMAL
25        berni(CELL) = pk(CELL)
26        ! from now on pk contains the vertically-averaged geopotential
27        pk(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
28      END_BLOCK
29    END_BLOCK
30  CASE(thermo_theta)
31    FORALL_CELLS()
32      ON_PRIMAL
33        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
34      END_BLOCK
35    END_BLOCK
36  CASE(thermo_entropy)
37    FORALL_CELLS()
38      ON_PRIMAL
39        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
40        berni(CELL) = berni(CELL) + pk(CELL)*(cpp-theta(CELL,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
41      END_BLOCK
42    END_BLOCK
43  CASE(thermo_variable_Cp)
44    ! thermodynamics with variable Cp
45    !           Cp(T) = Cp0 * (T/T0)^nu
46    ! =>            h = Cp(T).T/(nu+1)
47
48    FORALL_CELLS()
49      ON_PRIMAL
50        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
51        cp_ik = cpp*(pk(CELL)/Treff)**nu
52        berni(CELL) = berni(CELL) + pk(CELL)*(cp_ik/(nu+1.)-theta(CELL,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
53      END_BLOCK
54    END_BLOCK
55  CASE DEFAULT
56    PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo  ! FIXME
57    STOP
58  END SELECT
59!
60  FORALL_CELLS()
61    ON_EDGES
62      due = .5*(theta(CELL1,1)+theta(CELL2,1))*(pk(CELL2)-pk(CELL1)) + berni(CELL2)-berni(CELL1)
63      du(EDGE) = du(EDGE) - SIGN*due
64      u(EDGE) = u(EDGE) + tau*du(EDGE)
65    END_BLOCK
66  END_BLOCK 
67!
68END_BLOCK
69
70#endif END_DYSL
71
72  SUBROUTINE compute_caldyn_fast_unst(tau,theta,geopot, pk,berni,du,u)
73    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
74    USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
75         id_fast, primal_num, dual_num, edge_num, &
76         dual_deg, dual_edge, dual_ne, dual_vertex, &
77         up, down, left, right, Av, fv, Riv2
78    NUM, INTENT(IN) :: tau
79    FIELD_MASS   :: pk,berni  ! INOUT, OUT
80    FIELD_THETA  :: theta     ! IN
81    FIELD_GEOPOT :: geopot    ! IN
82    FIELD_U      :: u,du      ! INOUT,INOUT
83    DECLARE_INDICES
84    DECLARE_EDGES
85    NUM          :: due, cp_ik
86    START_TRACE(id_fast, 4,0,2) ! primal, dual, edge
87#include "../kernels_unst/caldyn_fast.k90"
88    STOP_TRACE
89  END SUBROUTINE compute_caldyn_fast_unst
90
91  SUBROUTINE compute_caldyn_fast_hex(tau,theta,geopot, pk,berni,du,u)
92    USE icosa
93    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
94    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
95    REAL(rstd),INTENT(IN)    :: geopot(iim*jjm,llm+1)
96    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
97    REAL(rstd),INTENT(INOUT) :: berni(iim*jjm,llm)  ! partial Bernoulli function
98    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
99    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! INOUT if tau>0
100    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
101
102    INTEGER :: ij,l
103    REAL(rstd) :: cp_ik, qv, temp, chi, log_p_preff, due, due_right, due_lup, due_ldown
104
105    CALL trace_start("compute_caldyn_fast")
106#include "../kernels_hex/caldyn_fast.k90"
107    CALL trace_end("compute_caldyn_fast")
108  END SUBROUTINE compute_caldyn_fast_hex
109
110  SUBROUTINE compute_caldyn_fast_manual(tau,theta,geopot, pk,berni,du,u)
111    USE icosa
112    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
113    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
114    REAL(rstd),INTENT(IN)    :: geopot(iim*jjm,llm+1)
115    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
116    REAL(rstd),INTENT(INOUT) :: berni(iim*jjm,llm)  ! partial Bernoulli function
117    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
118    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! INOUT if tau>0
119    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
120
121    INTEGER :: ij,l
122    REAL(rstd) :: cp_ik, qv, temp, chi, log_p_preff, due, due_right, due_lup, due_ldown
123
124    CALL trace_start("compute_caldyn_fast")
125
126    ! Compute Bernoulli term
127    IF(boussinesq) THEN
128       DO l=ll_begin,ll_end
129          !DIR$ SIMD
130          DO ij=ij_begin,ij_end         
131             berni(ij,l) = pk(ij,l)
132             ! from now on pk contains the vertically-averaged geopotential
133             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
134          END DO
135       END DO
136    ELSE ! compressible
137
138       DO l=ll_begin,ll_end
139          SELECT CASE(caldyn_thermo)
140          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
141             !DIR$ SIMD
142             DO ij=ij_begin,ij_end         
143                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
144             END DO
145          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)
146             !DIR$ SIMD
147             DO ij=ij_begin,ij_end         
148                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
149                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
150             END DO
151          CASE(thermo_moist) 
152             !DIR$ SIMD
153             DO ij=ij_begin,ij_end         
154                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
155                ! Bd = Phi + gibbs_d
156                ! Bv = Phi + gibbs_v
157                ! pk=temperature, theta=entropy
158                qv = theta(ij,l,2)
159                temp = pk(ij,l)
160                chi = log(temp/Treff)
161                log_p_preff = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
162                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
163                     + temp*(cpp*(1.-chi)+Rd*log_p_preff)
164                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
165                     + temp*(cppv*(1.-chi)+Rv*log_p_preff)
166            END DO
167          END SELECT
168       END DO
169
170    END IF ! Boussinesq/compressible
171
172!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
173    DO l=ll_begin,ll_end
174       IF(caldyn_thermo == thermo_moist) THEN
175          !DIR$ SIMD
176          DO ij=ij_begin,ij_end         
177             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
178                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
179                       *(pk(ij+t_right,l)-pk(ij,l))             &
180                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
181                       *(berniv(ij+t_right,l)-berniv(ij,l))
182             
183             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
184                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
185                       *(pk(ij+t_lup,l)-pk(ij,l))               &
186                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
187                       *(berniv(ij+t_lup,l)-berniv(ij,l))
188             
189             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
190                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
191                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
192                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
193                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
194             
195             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
196             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
197             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
198             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
199             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
200             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
201          END DO
202       ELSE
203          !DIR$ SIMD
204          DO ij=ij_begin,ij_end         
205             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
206                  *(pk(ij+t_right,l)-pk(ij,l))        &
207                  +  berni(ij+t_right,l)-berni(ij,l)
208             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
209                  *(pk(ij+t_lup,l)-pk(ij,l))          &
210                  +  berni(ij+t_lup,l)-berni(ij,l)
211             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
212                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
213                  +   berni(ij+t_ldown,l)-berni(ij,l)
214             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
215             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
216             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
217             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
218             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
219             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
220          END DO
221       END IF
222    END DO
223
224    CALL trace_end("compute_caldyn_fast")
225
226  END SUBROUTINE compute_caldyn_fast_manual
227
228END MODULE compute_caldyn_fast_mod
Note: See TracBrowser for help on using the repository browser.