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

Last change on this file since 1027 was 1027, checked in by dubos, 4 years ago

devel : towards conformity to F2008 standard
File size: 7.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!-------------- Wrappers for F2008 conformity -----------------
56
57  SUBROUTINE compute_caldyn_slow_hydro_hex(zero, u,rhodz,hv,Kv, berni, hflux,du)
58    LOGICAL, INTENT(IN)      :: zero
59    REAL(rstd),INTENT(IN)    :: u(:,:), rhodz(:,:), hv(:,:), Kv(:,:)
60    REAL(rstd), INTENT(OUT)  :: berni(:,:), hflux(:,:)
61    REAL(rstd),INTENT(INOUT) :: du(:,:)
62    CALL compute_caldyn_slow_hydro_hex_(zero, u,rhodz,hv,Kv, berni, hflux,du)
63  END SUBROUTINE compute_caldyn_slow_hydro_hex
64
65  SUBROUTINE compute_caldyn_slow_hydro_unst(zero, u,rhodz,hv,Kv, berni, hflux,du)
66    LOGICAL, INTENT(IN)      :: zero
67    REAL(rstd),INTENT(IN)    :: u(:,:), rhodz(:,:), hv(:,:), Kv(:,:)
68    REAL(rstd), INTENT(OUT)  :: berni(:,:), hflux(:,:)
69    REAL(rstd),INTENT(INOUT) :: du(:,:)
70    CALL compute_caldyn_slow_hydro_unst_(zero, u,rhodz,hv,Kv, berni, hflux,du)
71  END SUBROUTINE compute_caldyn_slow_hydro_unst
72
73!--------------------------------------------------------------
74
75  SUBROUTINE compute_caldyn_slow_hydro_unst_(zero, u,rhodz,hv,Kv, berni, hflux,du)
76    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
77    USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
78         id_slow_hydro, left, right, primal_deg, primal_edge
79    LOGICAL, INTENT(IN) :: zero
80    FIELD_MASS  :: rhodz, hv, Kv, berni  ! IN, IN, IN, BUF
81    FIELD_U     :: u,hflux,du ! IN, OUT, OUT
82    DECLARE_INDICES
83    DECLARE_EDGES
84    NUM :: ke, uu
85    START_TRACE(id_slow_hydro, 3,0,3)
86#include "../kernels_unst/caldyn_slow_hydro.k90"
87    STOP_TRACE
88  END SUBROUTINE compute_caldyn_slow_hydro_unst_
89
90  SUBROUTINE compute_caldyn_slow_hydro_hex_(zero, u,rhodz,hv,Kv, berni, hflux,du)
91    USE icosa
92    USE caldyn_vars_mod
93    LOGICAL, INTENT(IN) :: zero
94    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
95    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
96    REAL(rstd),INTENT(IN)  :: hv(2*iim*jjm,llm)   ! height/mass averaged to vertices
97    REAL(rstd),INTENT(IN)  :: Kv(2*iim*jjm,llm)   ! kinetic energy at vertices
98    REAL(rstd), INTENT(OUT) :: berni(iim*jjm,llm)  ! Bernoulli function
99    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
100    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
101    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
102    REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
103    INTEGER :: ij,l
104
105    CALL trace_start("compute_caldyn_slow_hydro")
106#define BERNI(ij,l) berni(ij,l)
107#include "../kernels_hex/caldyn_slow_hydro.k90"
108#undef BERNI
109    CALL trace_end("compute_caldyn_slow_hydro")   
110
111  END SUBROUTINE compute_caldyn_slow_hydro_hex_
112
113  SUBROUTINE compute_caldyn_slow_hydro_manual(zero, u,rhodz,hv,Kv, berni, hflux,du)
114    USE icosa
115    USE caldyn_vars_mod
116    LOGICAL, INTENT(IN) :: zero
117    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
118    REAL(rstd),INTENT(IN)  :: Kv(2*iim*jjm,llm)   ! kinetic energy at vertices
119    REAL(rstd),INTENT(IN)  :: hv(2*iim*jjm,llm)   ! height/mass averaged to vertices
120    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
121    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
122    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
123   
124    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
125    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
126    REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
127    INTEGER :: ij,l
128
129    CALL trace_start("compute_caldyn_slow_hydro")
130
131#define BERNI(ij) berni1(ij)
132
133    DO l = ll_begin, ll_end
134       !  Compute mass fluxes
135       IF (caldyn_conserv==conserv_energy) CALL test_message(req_qu) 
136
137       IF(caldyn_kinetic==kinetic_trisk) THEN
138          !DIR$ SIMD
139          DO ij=ij_begin_ext,ij_end_ext
140             uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
141             uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
142             uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
143             uu_right= uu_right*le_de(ij+u_right)
144             uu_lup  = uu_lup  *le_de(ij+u_lup)
145             uu_ldown= uu_ldown*le_de(ij+u_ldown)
146             hflux(ij+u_right,l)=uu_right
147             hflux(ij+u_lup,l)  =uu_lup
148             hflux(ij+u_ldown,l)=uu_ldown
149          ENDDO
150       ELSE ! mass flux deriving from consistent kinetic energy
151          !DIR$ SIMD
152          DO ij=ij_begin_ext,ij_end_ext
153             uu_right=0.5*(hv(ij+z_rup,l)+hv(ij+z_rdown,l))*u(ij+u_right,l)
154             uu_lup=0.5*(hv(ij+z_up,l)+hv(ij+z_lup,l))*u(ij+u_lup,l)
155             uu_ldown=0.5*(hv(ij+z_ldown,l)+hv(ij+z_down,l))*u(ij+u_ldown,l)
156             uu_right= uu_right*le_de(ij+u_right)
157             uu_lup  = uu_lup  *le_de(ij+u_lup)
158             uu_ldown= uu_ldown*le_de(ij+u_ldown)
159             hflux(ij+u_right,l)=uu_right
160             hflux(ij+u_lup,l)  =uu_lup
161             hflux(ij+u_ldown,l)=uu_ldown
162          ENDDO
163       END IF
164
165       ! Compute Bernoulli=kinetic energy
166       IF(caldyn_kinetic==kinetic_trisk) THEN
167          !DIR$ SIMD
168          DO ij=ij_begin,ij_end         
169             BERNI(ij) = &
170                  1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
171                                le_de(ij+u_rup)*u(ij+u_rup,l)**2     +    &
172                                le_de(ij+u_lup)*u(ij+u_lup,l)**2     +    &
173                                le_de(ij+u_left)*u(ij+u_left,l)**2   +    &
174                                le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
175                                le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
176          ENDDO
177       ELSE
178          !DIR$ SIMD
179          DO ij=ij_begin,ij_end
180             BERNI(ij) = Riv(ij,vup)   *Kv(ij+z_up,l)    + &
181                         Riv(ij,vlup)  *Kv(ij+z_lup,l)   + &
182                         Riv(ij,vldown)*Kv(ij+z_ldown,l) + &
183                         Riv(ij,vdown) *Kv(ij+z_down,l)  + &
184                         Riv(ij,vrdown)*Kv(ij+z_rdown,l) + &
185                         Riv(ij,vrup)  *Kv(ij+z_rup,l)
186          END DO
187       END IF
188       ! Compute du=-grad(Bernoulli)
189       IF(zero) THEN
190          !DIR$ SIMD
191          DO ij=ij_begin,ij_end
192             du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
193             du(ij+u_lup,l)   = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
194             du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
195          END DO
196       ELSE
197          !DIR$ SIMD
198          DO ij=ij_begin,ij_end
199             du(ij+u_right,l) = du(ij+u_right,l) + &
200                  ne_right*(BERNI(ij)-BERNI(ij+t_right))
201             du(ij+u_lup,l)   = du(ij+u_lup,l) + &
202                  ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
203             du(ij+u_ldown,l) = du(ij+u_ldown,l) + &
204                  ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
205          END DO
206       END IF
207    END DO
208
209#undef BERNI
210
211    CALL trace_end("compute_caldyn_slow_hydro")   
212  END SUBROUTINE compute_caldyn_slow_hydro_manual
213
214END MODULE compute_caldyn_slow_hydro_mod
Note: See TracBrowser for help on using the repository browser.