source: codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90 @ 735

Last change on this file since 735 was 735, checked in by dubos, 6 years ago

devel : Gassmann (2018) modification of TRiSK Coriolis discretization
File size: 6.5 KB
Line 
1MODULE caldyn_hevi_mod
2  USE icosa
3  USE transfert_mod
4  USE caldyn_vars_mod
5  USE caldyn_kernels_hevi_mod
6  USE caldyn_kernels_base_mod
7  IMPLICIT NONE
8  PRIVATE
9  PUBLIC caldyn_hevi
10
11CONTAINS
12 
13  SUBROUTINE caldyn_hevi(write_out,tau, f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
14       f_W, f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, &
15       f_du_slow, f_du_fast, f_dPhi_slow, f_dPhi_fast, f_dW_slow, f_dW_fast) 
16    USE icosa
17    USE observable_mod
18    USE disvert_mod, ONLY : caldyn_eta, eta_mass
19    USE vorticity_mod
20    USE kinetic_mod
21    USE theta2theta_rhodz_mod
22    USE wind_mod
23    USE mpipara
24    USE trace
25    USE omp_para
26    USE output_field_mod
27    USE checksum_mod
28    IMPLICIT NONE
29    LOGICAL,INTENT(IN)    :: write_out
30    REAL(rstd), INTENT(IN) :: tau
31    TYPE(t_field),POINTER :: f_phis(:)
32    TYPE(t_field),POINTER :: f_ps(:)
33    TYPE(t_field),POINTER :: f_mass(:)
34    TYPE(t_field),POINTER :: f_theta_rhodz(:)
35    TYPE(t_field),POINTER :: f_u(:)
36    TYPE(t_field),POINTER :: f_q(:)
37    TYPE(t_field),POINTER :: f_W(:)
38    TYPE(t_field),POINTER :: f_geopot(:)
39    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
40    TYPE(t_field) :: f_dps(:)
41    TYPE(t_field) :: f_dmass(:)
42    TYPE(t_field) :: f_dtheta_rhodz(:)
43    TYPE(t_field) :: f_du_slow(:)
44    TYPE(t_field) :: f_du_fast(:)
45    TYPE(t_field) :: f_dW_slow(:)
46    TYPE(t_field) :: f_dW_fast(:)
47    TYPE(t_field) :: f_dPhi_slow(:)
48    TYPE(t_field) :: f_dPhi_fast(:)
49   
50    REAL(rstd),POINTER :: ps(:), dps(:), phis(:)
51    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
52    REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
53    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:), Kv(:,:), hv(:,:)
54
55! temporary shared variable
56    REAL(rstd),POINTER  :: theta(:,:,:) 
57    REAL(rstd),POINTER  :: pk(:,:)
58    REAL(rstd),POINTER  :: geopot(:,:)
59    REAL(rstd),POINTER  :: convm(:,:) 
60    REAL(rstd),POINTER  :: wwuu(:,:)
61    REAL(rstd),POINTER  :: F_el(:,:), gradPhi2(:,:), w_il(:,:) , W_etadot(:,:), pres(:,:), m_il(:,:)
62    INTEGER :: ind
63    LOGICAL,SAVE :: first=.TRUE.
64!$OMP THREADPRIVATE(first)
65   
66    IF (first) THEN
67      first=.FALSE.
68      IF(caldyn_eta==eta_mass) THEN
69         CALL init_message(f_ps,req_i1,req_ps)
70      ELSE
71         CALL init_message(f_mass,req_i1,req_mass)
72      END IF
73      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
74      CALL init_message(f_u,req_e1_vect,req_u)
75      CALL init_message(f_qu,req_e1_scal,req_qu)
76      IF(caldyn_kinetic==kinetic_consistent) CALL init_message(f_Kv,req_z1_scal,req_Kv)
77      IF(.NOT.hydrostatic) THEN
78         CALL init_message(f_geopot,req_i1,req_geopot)
79         CALL init_message(f_w,req_i1,req_w)
80      END IF
81    ENDIF
82   
83    CALL trace_start("caldyn")
84   
85    IF(caldyn_eta==eta_mass) THEN
86       CALL send_message(f_ps,req_ps) ! COM00
87       CALL wait_message(req_ps) ! COM00
88    ELSE
89       CALL send_message(f_mass,req_mass) ! COM00
90       CALL wait_message(req_mass) ! COM00
91    END IF
92    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
93    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
94
95    IF(.NOT.hydrostatic) THEN
96       CALL send_message(f_geopot,req_geopot) ! COM03
97       CALL wait_message(req_geopot) ! COM03
98       CALL send_message(f_w,req_w) ! COM04
99       CALL wait_message(req_w) ! COM04
100    END IF
101   
102    DO ind=1,ndomain
103       IF (.NOT. assigned_domain(ind)) CYCLE
104       CALL swap_dimensions(ind)
105       CALL swap_geometry(ind)
106       ps=f_ps(ind)
107       theta_rhodz=f_theta_rhodz(ind)
108       mass=f_mass(ind)
109       theta = f_theta(ind)
110       CALL compute_theta(ps,theta_rhodz, mass,theta)
111       pk = f_pk(ind)
112       geopot = f_geopot(ind)
113       du=f_du_fast(ind)
114       IF(hydrostatic) THEN
115          du(:,:)=0.
116          CALL compute_geopot(mass,theta, ps,pk,geopot)
117       ELSE
118          phis = f_phis(ind)
119          W = f_W(ind)
120          dW = f_dW_fast(ind)
121          dPhi = f_dPhi_fast(ind)
122          ! reuse buffers
123          m_il = f_wil(ind)
124          pres = f_gradPhi2(ind)
125          CALL compute_caldyn_solver(tau,phis, mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
126       END IF
127       u=f_u(ind)
128       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot,du) ! computes du_fast and updates u
129    ENDDO
130   
131    CALL send_message(f_u,req_u) ! COM02
132    CALL wait_message(req_u)   ! COM02
133   
134    DO ind=1,ndomain
135       IF (.NOT. assigned_domain(ind)) CYCLE
136       CALL swap_dimensions(ind)
137       CALL swap_geometry(ind)
138       u=f_u(ind)
139       mass=f_mass(ind)
140       qu=f_qu(ind)
141       qv=f_qv(ind)
142       hv=f_hv(ind)
143       CALL compute_pvort_only(u,mass,qu,qv,hv)
144       IF(caldyn_kinetic==kinetic_consistent) THEN
145          Kv=f_Kv(ind)
146          CALL compute_caldyn_Kv(u,Kv)
147       END IF
148    ENDDO
149   
150    CALL send_message(f_qu,req_qu) ! COM03
151    CALL wait_message(req_qu) ! COM03
152
153    IF(caldyn_kinetic==kinetic_consistent) THEN
154       CALL transfert_request(f_hv,req_z1_scal)
155       CALL send_message(f_Kv,req_Kv)
156       CALL wait_message(req_Kv)
157    END IF
158
159    DO ind=1,ndomain
160       IF (.NOT. assigned_domain(ind)) CYCLE
161       CALL swap_dimensions(ind)
162       CALL swap_geometry(ind)
163       u=f_u(ind)
164       mass=f_mass(ind)
165       theta = f_theta(ind)
166       qu=f_qu(ind)
167       hflux=f_hflux(ind)
168       convm = f_dmass(ind)
169       dtheta_rhodz=f_dtheta_rhodz(ind)
170       du=f_du_slow(ind)
171
172       IF(hydrostatic) THEN
173          hv=f_hv(ind)
174          Kv=f_Kv(ind)
175          CALL compute_caldyn_slow_hydro(u,mass,hv, hflux,Kv,du, .TRUE.)
176       ELSE
177          W = f_W(ind)
178          dW = f_dW_slow(ind)
179          geopot = f_geopot(ind)
180          dPhi = f_dPhi_slow(ind)
181          F_el = f_Fel(ind)
182          gradPhi2 = f_gradPhi2(ind)
183          w_il = f_wil(ind)
184          CALL compute_caldyn_slow_NH(u,mass,geopot,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
185       END IF
186       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
187       IF(caldyn_eta==eta_mass) THEN
188          wflux=f_wflux(ind)
189          wwuu=f_wwuu(ind)
190          dps=f_dps(ind)
191          CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
192          IF(.NOT.hydrostatic) THEN
193             W_etadot=f_Wetadot(ind)
194             CALL compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
195          END IF
196       END IF
197    ENDDO
198   
199!$OMP BARRIER
200    !    CALL check_mass_conservation(f_ps,f_dps)
201    CALL trace_end("caldyn_hevi")
202!!$OMP BARRIER
203   
204  END SUBROUTINE caldyn_hevi
205
206END MODULE caldyn_hevi_mod
Note: See TracBrowser for help on using the repository browser.