source: codes/icosagcm/trunk/src/explicit_scheme.f90 @ 487

Last change on this file since 487 was 487, checked in by ymipsl, 8 years ago

Bad name for intel simd directives

YM
File size: 7.2 KB
Line 
1MODULE explicit_scheme_mod
2  USE euler_scheme_mod
3  USE prec
4  USE domain_mod
5  USE field_mod
6  IMPLICIT NONE
7  PRIVATE
8
9  REAL(rstd), DIMENSION(4), PARAMETER :: coef_rk4 = (/ .25, 1./3., .5, 1. /)
10  REAL(rstd), DIMENSION(5), PARAMETER :: coef_rk25 = (/ .25, 1./6., 3./8., .5, 1. /)
11
12  PUBLIC :: explicit_scheme
13CONTAINS
14
15  SUBROUTINE explicit_scheme(it, fluxt_zero)
16    USE time_mod
17    USE omp_para
18    USE caldyn_mod
19    USE trace
20    USE dimensions
21    USE geometry
22!    USE caldyn_gcm_mod, ONLY : req_ps, req_mass
23
24    REAL(rstd),POINTER :: q(:,:,:)
25    REAL(rstd),POINTER :: phis(:), ps(:) ,psm1(:), psm2(:), dps(:)
26    REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:)
27    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:)
28    REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:)
29    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
30    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
31    INTEGER :: it,stage
32    DO stage=1,nb_stage
33       !       CALL checksum(f_ps)
34       !       CALL checksum(f_theta_rhodz)
35       !       CALL checksum(f_mass)
36       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
37            f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, &
38            f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
39       !       CALL checksum(f_dps)
40       !       CALL checksum(f_dtheta_rhodz)
41       !       CALL checksum(f_dmass)
42       SELECT CASE (scheme)
43       CASE(euler)
44          CALL euler_scheme(.TRUE., fluxt_zero)
45       CASE (rk4)
46          CALL rk_scheme(stage, coef_rk4)
47       CASE (rk25)
48          CALL rk_scheme(stage, coef_rk25)
49       CASE (mlf)
50          CALL  leapfrog_matsuno_scheme(stage)
51       CASE DEFAULT
52          STOP
53       END SELECT
54    END DO
55
56  CONTAINS
57
58    SUBROUTINE RK_scheme(stage,coef)
59      USE disvert_mod
60      INTEGER :: ind, stage
61      REAL(rstd), INTENT(IN) :: coef(:)
62      REAL(rstd) :: tau
63      INTEGER :: i,j,ij,l
64
65      CALL trace_start("RK_scheme") 
66
67      tau = dt*coef(stage)
68
69      ! if mass coordinate, deal with ps first on one core
70      IF(caldyn_eta==eta_mass) THEN
71         IF (is_omp_first_level) THEN
72
73            DO ind=1,ndomain
74               IF (.NOT. assigned_domain(ind)) CYCLE
75               CALL swap_dimensions(ind)
76               CALL swap_geometry(ind)
77               ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind) 
78
79               IF (stage==1) THEN ! first stage : save model state in XXm1
80                  !DIR$ SIMD
81                  DO ij=ij_begin,ij_end
82                     psm1(ij)=ps(ij)
83                  ENDDO
84               ENDIF
85
86               ! updates are of the form x1 := x0 + tau*f(x1)
87               !DIR$ SIMD
88               DO ij=ij_begin,ij_end
89                  ps(ij)=psm1(ij)+tau*dps(ij)
90               ENDDO
91            ENDDO
92         ENDIF
93         !         CALL send_message(f_ps,req_ps)
94         !ym no overlap for now         
95         !         CALL wait_message(req_ps) 
96
97      ELSE ! Lagrangian coordinate, deal with mass
98         DO ind=1,ndomain
99            IF (.NOT. assigned_domain(ind)) CYCLE
100            CALL swap_dimensions(ind)
101            CALL swap_geometry(ind)
102            mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind)
103
104            IF (stage==1) THEN ! first stage : save model state in XXm1
105               DO l=ll_begin,ll_end
106                  !DIR$ SIMD
107                  DO ij=ij_begin,ij_end
108                     massm1(ij,l)=mass(ij,l)
109                  ENDDO
110               ENDDO
111            END IF
112
113            ! updates are of the form x1 := x0 + tau*f(x1)
114            DO l=ll_begin,ll_end
115               !DIR$ SIMD
116               DO ij=ij_begin,ij_end
117                  mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
118               ENDDO
119            ENDDO
120         END DO
121         !         CALL send_message(f_mass,req_mass)
122         !ym no overlap for now         
123         !         CALL wait_message(req_mass) 
124
125      END IF
126
127      ! now deal with other prognostic variables
128      DO ind=1,ndomain
129         IF (.NOT. assigned_domain(ind)) CYCLE
130         CALL swap_dimensions(ind)
131         CALL swap_geometry(ind)
132         u=f_u(ind)      ; du=f_du(ind)      ; um1=f_um1(ind) 
133         theta_rhodz=f_theta_rhodz(ind)
134         theta_rhodzm1=f_theta_rhodzm1(ind)
135         dtheta_rhodz=f_dtheta_rhodz(ind)
136
137         IF (stage==1) THEN ! first stage : save model state in XXm1
138            DO l=ll_begin,ll_end
139               !DIR$ SIMD
140               DO ij=ij_begin,ij_end
141                  um1(ij+u_right,l)=u(ij+u_right,l)
142                  um1(ij+u_lup,l)=u(ij+u_lup,l)
143                  um1(ij+u_ldown,l)=u(ij+u_ldown,l)
144                  theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
145               ENDDO
146            ENDDO
147         END IF
148
149         DO l=ll_begin,ll_end
150            !DIR$ SIMD
151            DO ij=ij_begin,ij_end
152               u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l)
153               u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
154               u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
155               theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
156            ENDDO
157         ENDDO
158
159         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
160            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
161            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
162            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
163         END IF
164      END DO
165
166      CALL trace_end("RK_scheme")
167
168    END SUBROUTINE RK_scheme
169
170    SUBROUTINE leapfrog_matsuno_scheme(stage)
171      IMPLICIT NONE
172      INTEGER :: ind, stage
173      REAL :: tau
174
175      CALL trace_start("leapfrog_matsuno_scheme")
176
177      tau = dt/nb_stage
178      DO ind=1,ndomain
179         IF (.NOT. assigned_domain(ind)) CYCLE
180         CALL swap_dimensions(ind)
181         CALL swap_geometry(ind)
182
183         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
184         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
185         psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
186         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
187
188
189         IF (stage==1) THEN
190            psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
191            ps(:)=psm1(:)+tau*dps(:)
192            u(:,:)=um1(:,:)+tau*du(:,:)
193            theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
194
195         ELSE IF (stage==2) THEN
196
197            ps(:)=psm1(:)+tau*dps(:)
198            u(:,:)=um1(:,:)+tau*du(:,:)
199            theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
200
201            psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
202            psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
203
204         ELSE
205
206            ps(:)=psm2(:)+2*tau*dps(:)
207            u(:,:)=um2(:,:)+2*tau*du(:,:)
208            theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
209
210            psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
211            psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
212
213         ENDIF
214
215      ENDDO
216      CALL trace_end("leapfrog_matsuno_scheme")
217
218    END SUBROUTINE leapfrog_matsuno_scheme
219
220  END SUBROUTINE Explicit_scheme
221
222END MODULE explicit_scheme_mod
Note: See TracBrowser for help on using the repository browser.