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

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

devel : interfaces for caldyn_fast and caldyn_slow_hydro
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
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, primal_num, dual_num, edge_num, &
74         dual_deg, dual_edge, dual_ne, dual_vertex, &
75         up, down, left, right, Av, fv, Riv2
76    NUM, INTENT(IN) :: tau
77    FIELD_MASS   :: pk,berni  ! INOUT, OUT
78    FIELD_THETA  :: theta     ! IN
79    FIELD_GEOPOT :: geopot    ! IN
80    FIELD_U      :: u,du      ! INOUT,INOUT
81    DECLARE_INDICES
82    DECLARE_EDGES
83    NUM          :: due, cp_ik, Phi_ik
84    START_TRACE(id_fast, 4,0,2) ! primal, dual, edge
85#include "../kernels_unst/caldyn_fast.k90"
86    STOP_TRACE
87  END SUBROUTINE compute_caldyn_fast_unst
88
89  SUBROUTINE compute_caldyn_fast_hex(tau,theta,geopot, pk,berni,du,u)
90    USE icosa
91    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
92    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
93    REAL(rstd),INTENT(IN)    :: geopot(iim*jjm,llm+1)
94    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
95    REAL(rstd),INTENT(INOUT) :: berni(iim*jjm,llm)  ! partial Bernoulli function
96    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
97    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! INOUT if tau>0
98    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
99
100    INTEGER :: ij,l
101    REAL(rstd) :: due, cp_ik, Phi_ik
102
103    CALL trace_start("compute_caldyn_fast")
104#include "../kernels_hex/caldyn_fast.k90"
105    CALL trace_end("compute_caldyn_fast")
106   
107  END SUBROUTINE compute_caldyn_fast_hex
108
109  SUBROUTINE compute_caldyn_fast_manual(tau,theta,geopot, pk,berni,du,u)
110    USE icosa
111    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
112    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
113    REAL(rstd),INTENT(IN)    :: geopot(iim*jjm,llm+1)
114    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
115    REAL(rstd),INTENT(INOUT) :: berni(iim*jjm,llm)  ! partial Bernoulli function
116    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
117    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! INOUT if tau>0
118    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
119
120    INTEGER :: ij,l
121    REAL(rstd) :: cp_ik, qv, temp, chi, log_p_preff, due, due_right, due_lup, due_ldown
122
123    CALL trace_start("compute_caldyn_fast")
124
125    ! Compute Bernoulli term
126    IF(boussinesq) THEN
127       DO l=ll_begin,ll_end
128          !DIR$ SIMD
129          DO ij=ij_begin,ij_end         
130             berni(ij,l) = pk(ij,l)
131             ! from now on pk contains the vertically-averaged geopotential
132             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
133          END DO
134       END DO
135    ELSE ! compressible
136
137       DO l=ll_begin,ll_end
138          SELECT CASE(caldyn_thermo)
139          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
140             !DIR$ SIMD
141             DO ij=ij_begin,ij_end         
142                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
143             END DO
144          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)
145             !DIR$ SIMD
146             DO ij=ij_begin,ij_end         
147                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
148                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
149             END DO
150          CASE(thermo_moist) 
151             !DIR$ SIMD
152             DO ij=ij_begin,ij_end         
153                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
154                ! Bd = Phi + gibbs_d
155                ! Bv = Phi + gibbs_v
156                ! pk=temperature, theta=entropy
157                qv = theta(ij,l,2)
158                temp = pk(ij,l)
159                chi = log(temp/Treff)
160                log_p_preff = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
161                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
162                     + temp*(cpp*(1.-chi)+Rd*log_p_preff)
163                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
164                     + temp*(cppv*(1.-chi)+Rv*log_p_preff)
165            END DO
166          END SELECT
167       END DO
168
169    END IF ! Boussinesq/compressible
170
171!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
172    DO l=ll_begin,ll_end
173       IF(caldyn_thermo == thermo_moist) THEN
174          !DIR$ SIMD
175          DO ij=ij_begin,ij_end         
176             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
177                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
178                       *(pk(ij+t_right,l)-pk(ij,l))             &
179                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
180                       *(berniv(ij+t_right,l)-berniv(ij,l))
181             
182             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
183                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
184                       *(pk(ij+t_lup,l)-pk(ij,l))               &
185                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
186                       *(berniv(ij+t_lup,l)-berniv(ij,l))
187             
188             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
189                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
190                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
191                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
192                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
193             
194             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
195             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
196             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
197             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
198             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
199             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
200          END DO
201       ELSE
202          !DIR$ SIMD
203          DO ij=ij_begin,ij_end         
204             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
205                  *(pk(ij+t_right,l)-pk(ij,l))        &
206                  +  berni(ij+t_right,l)-berni(ij,l)
207             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
208                  *(pk(ij+t_lup,l)-pk(ij,l))          &
209                  +  berni(ij+t_lup,l)-berni(ij,l)
210             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
211                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
212                  +   berni(ij+t_ldown,l)-berni(ij,l)
213             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
214             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
215             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
216             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
217             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
218             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
219          END DO
220       END IF
221    END DO
222
223    CALL trace_end("compute_caldyn_fast")
224
225  END SUBROUTINE compute_caldyn_fast_manual
226
227END MODULE compute_caldyn_fast_mod
Note: See TracBrowser for help on using the repository browser.