source: codes/icosagcm/trunk/src/time/euler_scheme.F90 @ 1057

Last change on this file since 1057 was 953, checked in by adurocher, 5 years ago

trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )
File size: 8.4 KB
Line 
1MODULE euler_scheme_mod
2  USE field_mod
3  USE abort_mod
4  IMPLICIT NONE
5  PRIVATE
6
7  TYPE(t_field),POINTER,SAVE,PUBLIC :: f_geopot(:), f_q(:), &
8       f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:), & ! current and previous time steps + tendency of mass,
9       f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:), &           ! column-integrated mass
10       f_u(:),f_um1(:),f_um2(:), f_du(:), f_W(:), &                   ! 'horizontal' and vertical (NH) momentum
11       f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:), & ! mass-weighted theta/entropy
12       f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:), &            ! accumulated mass fluxes
13       f_du_slow(:,:), f_du_fast(:,:), &                              ! slow/fast trend arrays
14       f_dmass_slow(:,:), f_dps_slow(:,:), f_dtheta_rhodz_slow(:,:), &! for HEVI time scheme
15       f_dPhi_slow(:,:), f_dPhi_fast(:,:), &                          ! geopotential tendencies (NH)
16       f_dW_slow(:,:), f_dW_fast(:,:)                                 ! vertical momentum tendencies (NH)
17
18  INTEGER, PARAMETER, PUBLIC :: explicit=1, hevi=2, euler=1, rk4=2, mlf=3, rk25=4, ark23=6, ark33=7
19
20  INTEGER,SAVE, PUBLIC :: nb_stage, matsuno_period, scheme, scheme_family
21  !$OMP THREADPRIVATE(nb_stage, matsuno_period, scheme, scheme_family)
22
23  PUBLIC :: euler_scheme, accumulate_fluxes, legacy_to_DEC, DEC_to_legacy
24CONTAINS
25
26  SUBROUTINE Euler_scheme(with_dps,fluxt_zero)
27    USE icosa
28    USE disvert_mod
29    USE omp_para
30    USE trace
31    LOGICAL :: with_dps
32    LOGICAL, OPTIONAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
33    REAL(rstd),POINTER :: ps(:), dps(:)
34    REAL(rstd),POINTER :: u(:,:) , du(:,:)
35    REAL(rstd),POINTER :: mass(:,:), dmass(:,:)
36    REAL(rstd),POINTER :: theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
37    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
38    INTEGER :: ind
39    INTEGER :: ij,l
40    CALL trace_start("Euler_scheme") 
41
42    DO ind=1,ndomain
43       IF (.NOT. assigned_domain(ind)) CYCLE
44       CALL swap_dimensions(ind)
45       CALL swap_geometry(ind)
46
47       IF(with_dps) THEN ! update ps/mass
48          CALL abort_acc("Euler_scheme/with_dps")
49          IF(caldyn_eta==eta_mass) THEN ! update ps
50             ps=f_ps(ind) ; dps=f_dps(ind) ;             
51             IF (is_omp_first_level) THEN
52                !DIR$ SIMD
53                DO ij=ij_begin,ij_end
54                   ps(ij)=ps(ij)+dt*dps(ij)
55                ENDDO
56             ENDIF
57          ELSE ! update mass
58             mass=f_mass(ind) ; dmass=f_dmass(ind) ;             
59             DO l=ll_begin,ll_end
60                !DIR$ SIMD
61                DO ij=ij_begin,ij_end
62                   mass(ij,l)=mass(ij,l)+dt*dmass(ij,l)
63                ENDDO
64             END DO
65          END IF
66
67          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
68          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
69          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
70       END IF ! update ps/mass
71
72       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
73       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
74
75       CALL compute_euler_scheme(u, du, theta_rhodz, dtheta_rhodz)
76    ENDDO
77
78    CALL trace_end("Euler_scheme")
79
80    CONTAINS
81
82    SUBROUTINE compute_euler_scheme(u, du, theta_rhodz, dtheta_rhodz)
83       REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)
84       REAL(rstd),INTENT(IN)    :: du(iim*3*jjm,llm)
85       REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm,nqdyn)
86       REAL(rstd),INTENT(IN)    :: dtheta_rhodz(iim*jjm,llm,nqdyn)
87
88       !$acc data present(theta_rhodz(:,:,:), u(:,:),du(:,:), dtheta_rhodz(:,:,:)) async
89
90       !$acc parallel loop async
91       DO l=ll_begin,ll_end
92          !$acc loop
93          !DIR$ SIMD
94          DO ij=ij_begin,ij_end
95             u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l)
96             u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
97             u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
98             theta_rhodz(ij,l,1)=theta_rhodz(ij,l,1)+dt*dtheta_rhodz(ij,l,1)
99          ENDDO
100       ENDDO
101       
102       !$acc end data
103    END SUBROUTINE compute_euler_scheme
104
105  END SUBROUTINE Euler_scheme
106
107  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
108    USE dimensions
109    USE omp_para
110    USE disvert_mod
111    IMPLICIT NONE
112    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
113    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
114    REAL(rstd), INTENT(IN) :: tau
115    LOGICAL, INTENT(INOUT) :: fluxt_zero
116    INTEGER :: l,ij
117
118    !$acc data present(hflux(:,:), wflux(:,:), hfluxt(:,:), wfluxt(:,:)) async
119
120    IF(fluxt_zero) THEN
121
122       fluxt_zero=.FALSE.
123       !$acc parallel loop async
124       DO l=ll_begin,ll_end
125          !$acc loop
126          !DIR$ SIMD
127          DO ij=ij_begin_ext,ij_end_ext
128             hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l)
129             hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l)
130             hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l)
131          ENDDO
132       ENDDO
133
134       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
135          !$acc parallel loop async
136          DO l=ll_begin,ll_endp1
137             !$acc loop
138             !DIR$ SIMD
139             DO ij=ij_begin,ij_end
140                wfluxt(ij,l) = tau*wflux(ij,l)
141             ENDDO
142          ENDDO
143       END IF
144
145    ELSE
146       !$acc parallel loop async
147       DO l=ll_begin,ll_end
148          !$acc loop
149          !DIR$ SIMD
150          DO ij=ij_begin_ext,ij_end_ext
151             hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l)
152             hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l)
153             hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l)
154          ENDDO
155       ENDDO
156
157       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
158          !$acc parallel loop async
159          DO l=ll_begin,ll_endp1
160             !$acc loop
161             !DIR$ SIMD
162             DO ij=ij_begin,ij_end
163                wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
164             ENDDO
165          ENDDO
166       END IF
167    END IF
168    !$acc end data
169  END SUBROUTINE accumulate_fluxes
170
171  SUBROUTINE legacy_to_DEC(f_ps, f_u)
172    USE icosa
173    USE disvert_mod
174    USE omp_para
175    USE trace
176    TYPE(t_field),POINTER :: f_ps(:), f_u(:)
177    REAL(rstd), POINTER :: ps(:), u(:,:)
178    INTEGER :: ind,ij,l
179
180    CALL trace_start("legacy_to_DEC")
181
182    DO ind=1,ndomain
183       IF (.NOT. assigned_domain(ind)) CYCLE
184       CALL swap_dimensions(ind)
185       CALL swap_geometry(ind)
186
187       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN ! update ps
188          ps=f_ps(ind)
189          !$acc parallel loop async default(present)
190          !DIR$ SIMD
191          DO ij=ij_begin,ij_end
192             ps(ij)=(ps(ij)-ptop)/g ! convert ps to column-integrated mass
193          ENDDO
194       END IF
195
196       u=f_u(ind)
197       !$acc parallel loop async default(present)
198       DO l=ll_begin,ll_end
199          !$acc loop
200          !DIR$ SIMD
201          DO ij=ij_begin,ij_end
202             u(ij+u_right,l)=u(ij+u_right,l)*de(ij+u_right)
203             u(ij+u_lup,l)=u(ij+u_lup,l)*de(ij+u_lup)
204             u(ij+u_ldown,l)=u(ij+u_ldown,l)*de(ij+u_ldown)
205          ENDDO
206       ENDDO
207    ENDDO
208
209    CALL trace_end("legacy_to_DEC") 
210
211  END SUBROUTINE Legacy_to_DEC
212
213  SUBROUTINE DEC_to_legacy(f_ps, f_u)
214    USE icosa
215    USE disvert_mod
216    USE omp_para
217    USE trace
218    TYPE(t_field),POINTER :: f_ps(:), f_u(:)
219    REAL(rstd), POINTER :: ps(:), u(:,:)
220    INTEGER :: ind,ij,l
221    CALL trace_start("legacy_to_DEC")
222
223    DO ind=1,ndomain
224       IF (.NOT. assigned_domain(ind)) CYCLE
225       CALL swap_dimensions(ind)
226       CALL swap_geometry(ind)
227
228       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN
229          ps=f_ps(ind)
230          !$acc parallel loop async default(present)
231          !DIR$ SIMD
232          DO ij=ij_begin,ij_end
233             ps(ij)=ptop+ps(ij)*g ! convert column-integrated mass to ps
234          ENDDO
235       ENDIF
236       
237       u=f_u(ind)
238       !$acc parallel loop async default(present)
239       DO l=ll_begin,ll_end
240          !$acc loop
241          !DIR$ SIMD
242          DO ij=ij_begin,ij_end
243             u(ij+u_right,l)=u(ij+u_right,l)/de(ij+u_right)
244             u(ij+u_lup,l)=u(ij+u_lup,l)/de(ij+u_lup)
245             u(ij+u_ldown,l)=u(ij+u_ldown,l)/de(ij+u_ldown)
246          ENDDO
247       ENDDO
248    ENDDO
249
250    CALL trace_end("DEC_to_legacy") 
251  END SUBROUTINE DEC_to_legacy
252
253END MODULE euler_scheme_mod
Note: See TracBrowser for help on using the repository browser.