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

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

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