source: codes/icosagcm/trunk/src/caldyn_hevi.f90 @ 361

Last change on this file since 361 was 361, checked in by dubos, 9 years ago

ARK2.3 HEVI time stepping - tested with DCMIP4 only

File size: 4.4 KB
Line 
1MODULE caldyn_hevi_mod
2  USE icosa
3  USE transfert_mod
4  USE caldyn_kernels_mod
5  USE caldyn_gcm_mod
6
7  IMPLICIT NONE
8  PRIVATE
9
10  PUBLIC caldyn_hevi
11
12CONTAINS
13 
14  SUBROUTINE caldyn_hevi(write_out,tau, f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
15       f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du_slow, f_du_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_geopot(:)
38    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
39    TYPE(t_field) :: f_dps(:)
40    TYPE(t_field) :: f_dmass(:)
41    TYPE(t_field) :: f_dtheta_rhodz(:)
42    TYPE(t_field) :: f_du_slow(:)
43    TYPE(t_field) :: f_du_fast(:)
44   
45    REAL(rstd),POINTER :: ps(:), dps(:)
46    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
47    REAL(rstd),POINTER :: du(:,:), hflux(:,:), wflux(:,:)
48    REAL(rstd),POINTER :: u(:,:), qu(:,:), qv(:,:)
49
50! temporary shared variable
51    REAL(rstd),POINTER  :: theta(:,:) 
52    REAL(rstd),POINTER  :: pk(:,:)
53    REAL(rstd),POINTER  :: geopot(:,:)
54    REAL(rstd),POINTER  :: convm(:,:) 
55    REAL(rstd),POINTER  :: wwuu(:,:)
56       
57    INTEGER :: ind
58    LOGICAL,SAVE :: first=.TRUE.
59!$OMP THREADPRIVATE(first)
60   
61    IF (first) THEN
62      first=.FALSE.
63      IF(caldyn_eta==eta_mass) THEN
64         CALL init_message(f_ps,req_i1,req_ps)
65      ELSE
66         CALL init_message(f_mass,req_i1,req_mass)
67      END IF
68      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
69      CALL init_message(f_u,req_e1_vect,req_u)
70      CALL init_message(f_qu,req_e1_scal,req_qu)
71    ENDIF
72   
73    CALL trace_start("caldyn")
74   
75    IF(caldyn_eta==eta_mass) THEN
76       CALL send_message(f_ps,req_ps) ! COM00
77    ELSE
78       CALL send_message(f_mass,req_mass) ! COM00
79    END IF
80
81    ! Block below moved from caldyn_pvort when splitting it into caldyn_theta and caldyn_pvort_only
82    IF(caldyn_eta==eta_mass) THEN
83       CALL wait_message(req_ps) ! COM00
84    ELSE
85       CALL wait_message(req_mass) ! COM00
86    END IF
87
88    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
89    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
90   
91    DO ind=1,ndomain
92       IF (.NOT. assigned_domain(ind)) CYCLE
93       CALL swap_dimensions(ind)
94       CALL swap_geometry(ind)
95       ps=f_ps(ind)
96       u=f_u(ind)
97       du=f_du_fast(ind)
98       theta_rhodz=f_theta_rhodz(ind)
99       mass=f_mass(ind)
100       theta = f_theta(ind)
101       pk = f_pk(ind)
102       geopot = f_geopot(ind) 
103       CALL compute_theta(ps,theta_rhodz, mass,theta)
104       CALL compute_geopot(ps,mass,theta, pk,geopot)
105       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot, du) ! computes du_fast and updates u
106    ENDDO
107   
108    CALL send_message(f_u,req_u) ! COM02
109    CALL wait_message(req_u)   ! COM02
110   
111    DO ind=1,ndomain
112       IF (.NOT. assigned_domain(ind)) CYCLE
113       CALL swap_dimensions(ind)
114       CALL swap_geometry(ind)
115       u=f_u(ind)
116       mass=f_mass(ind)
117       qu=f_qu(ind)
118       qv=f_qv(ind)
119       CALL compute_pvort_only(u,mass,qu,qv)
120    ENDDO
121   
122    CALL send_message(f_qu,req_qu) ! COM03
123   
124    DO ind=1,ndomain
125       IF (.NOT. assigned_domain(ind)) CYCLE
126       CALL swap_dimensions(ind)
127       CALL swap_geometry(ind)
128       u=f_u(ind)
129       mass=f_mass(ind)
130       theta = f_theta(ind)
131       qu=f_qu(ind)
132       hflux=f_hflux(ind)
133       convm = f_dmass(ind)
134       dtheta_rhodz=f_dtheta_rhodz(ind)
135       du=f_du_slow(ind)
136       CALL compute_caldyn_slow(u,mass,qu,theta, hflux,convm,dtheta_rhodz,du) ! COM03
137       IF(caldyn_eta==eta_mass) THEN
138          wflux=f_wflux(ind)
139          wwuu=f_wwuu(ind)
140          dps=f_dps(ind)
141          CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
142       END IF
143    ENDDO
144   
145!$OMP BARRIER
146    !    CALL check_mass_conservation(f_ps,f_dps)
147    CALL trace_end("caldyn_hevi")
148!!$OMP BARRIER
149   
150  END SUBROUTINE caldyn_hevi
151
152END MODULE caldyn_hevi_mod
Note: See TracBrowser for help on using the repository browser.