source: codes/icosagcm/trunk/src/euler_scheme.f90 @ 364

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

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