source: codes/icosagcm/trunk/src/time/explicit_scheme.f90 @ 1015

Last change on this file since 1015 was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

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