source: codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_hydro.F90 @ 939

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

devel : cleanup USE data_unstructured_mod
File size: 6.7 KB
Line 
1MODULE compute_caldyn_slow_hydro_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_slow_hydro_unst, &
14       compute_caldyn_slow_hydro_hex, compute_caldyn_slow_hydro_manual
15
16CONTAINS
17
18#ifdef BEGIN_DYSL
19
20KERNEL(caldyn_slow_hydro)
21  FORALL_CELLS_EXT()
22    ON_EDGES
23      uu = .5*(rhodz(CELL1)+rhodz(CELL2))*u(EDGE)
24      hflux(EDGE) = uu*LE_DE
25    END_BLOCK
26  END_BLOCK
27
28  FORALL_CELLS()
29    ON_PRIMAL
30      ke=0.d0
31      FORALL_EDGES
32        ke = ke + LE_DE*u(EDGE)**2
33      END_BLOCK
34      BERNI(CELL)=ke*(.25/AI)
35    END_BLOCK
36  END_BLOCK
37  IF(zero) THEN
38    FORALL_CELLS()
39      ON_EDGES
40        du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2)) ! minus gradient
41      END_BLOCK
42    END_BLOCK
43  ELSE
44    FORALL_CELLS()
45      ON_EDGES
46        du(EDGE) = du(EDGE) + SIGN*(berni(CELL1)-berni(CELL2)) ! minus gradient
47      END_BLOCK
48    END_BLOCK
49  END IF
50
51END_BLOCK
52
53#endif END_DYSL
54
55  SUBROUTINE compute_caldyn_slow_hydro_unst(zero, u,rhodz,hv,Kv, berni, hflux,du)
56    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
57    USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
58         id_slow_hydro, left, right, primal_deg, primal_edge
59    LOGICAL, INTENT(IN) :: zero
60    FIELD_MASS  :: rhodz, hv, Kv, berni  ! IN, IN, IN, BUF
61    FIELD_U     :: u,hflux,du ! IN, OUT, OUT
62    DECLARE_INDICES
63    DECLARE_EDGES
64    NUM :: ke, uu
65    START_TRACE(id_slow_hydro, 3,0,3)
66#include "../kernels_unst/caldyn_slow_hydro.k90"
67    STOP_TRACE
68  END SUBROUTINE compute_caldyn_slow_hydro_unst
69
70  SUBROUTINE compute_caldyn_slow_hydro_hex(zero, u,rhodz,hv,Kv, berni, hflux,du)
71    USE icosa
72    USE caldyn_vars_mod
73    LOGICAL, INTENT(IN) :: zero
74    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
75    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
76    REAL(rstd),INTENT(IN)  :: hv(2*iim*jjm,llm)   ! height/mass averaged to vertices
77    REAL(rstd),INTENT(IN)  :: Kv(2*iim*jjm,llm)   ! kinetic energy at vertices
78    REAL(rstd), INTENT(OUT) :: berni(iim*jjm,llm)  ! Bernoulli function
79    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
80    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
81    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
82    REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
83    INTEGER :: ij,l
84
85    CALL trace_start("compute_caldyn_slow_hydro")
86#define BERNI(ij,l) berni(ij,l)
87#include "../kernels_hex/caldyn_slow_hydro.k90"
88#undef BERNI
89    CALL trace_end("compute_caldyn_slow_hydro")   
90
91  END SUBROUTINE compute_caldyn_slow_hydro_hex
92
93  SUBROUTINE compute_caldyn_slow_hydro_manual(zero, u,rhodz,hv,Kv, berni, hflux,du)
94    USE icosa
95    USE caldyn_vars_mod
96    LOGICAL, INTENT(IN) :: zero
97    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
98    REAL(rstd),INTENT(IN)  :: Kv(2*iim*jjm,llm)   ! kinetic energy at vertices
99    REAL(rstd),INTENT(IN)  :: hv(2*iim*jjm,llm)   ! height/mass averaged to vertices
100    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
101    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
102    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
103   
104    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
105    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
106    REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
107    INTEGER :: ij,l
108
109    CALL trace_start("compute_caldyn_slow_hydro")
110
111#define BERNI(ij) berni1(ij)
112
113    DO l = ll_begin, ll_end
114       !  Compute mass fluxes
115       IF (caldyn_conserv==conserv_energy) CALL test_message(req_qu) 
116
117       IF(caldyn_kinetic==kinetic_trisk) THEN
118          !DIR$ SIMD
119          DO ij=ij_begin_ext,ij_end_ext
120             uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
121             uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
122             uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
123             uu_right= uu_right*le_de(ij+u_right)
124             uu_lup  = uu_lup  *le_de(ij+u_lup)
125             uu_ldown= uu_ldown*le_de(ij+u_ldown)
126             hflux(ij+u_right,l)=uu_right
127             hflux(ij+u_lup,l)  =uu_lup
128             hflux(ij+u_ldown,l)=uu_ldown
129          ENDDO
130       ELSE ! mass flux deriving from consistent kinetic energy
131          !DIR$ SIMD
132          DO ij=ij_begin_ext,ij_end_ext
133             uu_right=0.5*(hv(ij+z_rup,l)+hv(ij+z_rdown,l))*u(ij+u_right,l)
134             uu_lup=0.5*(hv(ij+z_up,l)+hv(ij+z_lup,l))*u(ij+u_lup,l)
135             uu_ldown=0.5*(hv(ij+z_ldown,l)+hv(ij+z_down,l))*u(ij+u_ldown,l)
136             uu_right= uu_right*le_de(ij+u_right)
137             uu_lup  = uu_lup  *le_de(ij+u_lup)
138             uu_ldown= uu_ldown*le_de(ij+u_ldown)
139             hflux(ij+u_right,l)=uu_right
140             hflux(ij+u_lup,l)  =uu_lup
141             hflux(ij+u_ldown,l)=uu_ldown
142          ENDDO
143       END IF
144
145       ! Compute Bernoulli=kinetic energy
146       IF(caldyn_kinetic==kinetic_trisk) THEN
147          !DIR$ SIMD
148          DO ij=ij_begin,ij_end         
149             BERNI(ij) = &
150                  1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
151                                le_de(ij+u_rup)*u(ij+u_rup,l)**2     +    &
152                                le_de(ij+u_lup)*u(ij+u_lup,l)**2     +    &
153                                le_de(ij+u_left)*u(ij+u_left,l)**2   +    &
154                                le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
155                                le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
156          ENDDO
157       ELSE
158          !DIR$ SIMD
159          DO ij=ij_begin,ij_end
160             BERNI(ij) = Riv(ij,vup)   *Kv(ij+z_up,l)    + &
161                         Riv(ij,vlup)  *Kv(ij+z_lup,l)   + &
162                         Riv(ij,vldown)*Kv(ij+z_ldown,l) + &
163                         Riv(ij,vdown) *Kv(ij+z_down,l)  + &
164                         Riv(ij,vrdown)*Kv(ij+z_rdown,l) + &
165                         Riv(ij,vrup)  *Kv(ij+z_rup,l)
166          END DO
167       END IF
168       ! Compute du=-grad(Bernoulli)
169       IF(zero) THEN
170          !DIR$ SIMD
171          DO ij=ij_begin,ij_end
172             du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
173             du(ij+u_lup,l)   = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
174             du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
175          END DO
176       ELSE
177          !DIR$ SIMD
178          DO ij=ij_begin,ij_end
179             du(ij+u_right,l) = du(ij+u_right,l) + &
180                  ne_right*(BERNI(ij)-BERNI(ij+t_right))
181             du(ij+u_lup,l)   = du(ij+u_lup,l) + &
182                  ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
183             du(ij+u_ldown,l) = du(ij+u_ldown,l) + &
184                  ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
185          END DO
186       END IF
187    END DO
188
189#undef BERNI
190
191    CALL trace_end("compute_caldyn_slow_hydro")   
192  END SUBROUTINE compute_caldyn_slow_hydro_manual
193
194END MODULE compute_caldyn_slow_hydro_mod
Note: See TracBrowser for help on using the repository browser.