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

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

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