source: codes/icosagcm/devel/src/time/explicit_scheme.f90 @ 714

Last change on this file since 714 was 533, checked in by dubos, 7 years ago

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