source: codes/icosagcm/devel/src/time/hevi_scheme.f90 @ 732

Last change on this file since 732 was 732, checked in by dubos, 6 years ago

devel : more cleanup and reorganization in dynamics/

File size: 5.4 KB
RevLine 
[360]1MODULE hevi_scheme_mod
2  USE prec
3  USE domain_mod
4  USE field_mod
[362]5  USE euler_scheme_mod
[732]6  USE caldyn_vars_mod
[360]7  IMPLICIT NONE
8  PRIVATE
9
10  REAL(rstd), SAVE :: wj(3), bjl(3,3), cjl(3,3), taujj(3)
11  REAL(rstd), PARAMETER, DIMENSION(3) :: zero = (/0.,0.,0./)
12
13  PUBLIC :: set_coefs_ark23, set_coefs_ark33, hevi_scheme
14
15CONTAINS
16
17  SUBROUTINE set_coefs_ark23(dt)
18    ! ARK2 scheme by Giraldo, Kelly, Constantinescu 2013
19    ! See Weller et al., 2013 - ARK2 scheme Fig. 2
20    REAL(rstd) :: dt
[371]21    REAL(rstd), PARAMETER :: delta=.5/SQRT(2.), gamma=1.-2.*delta
22!    REAL(rstd), PARAMETER :: alpha=(3.+SQRT(8.))/6. ! original value in GKC2013
23    REAL(rstd), PARAMETER :: alpha=0.7
[360]24    REAL(rstd), PARAMETER, DIMENSION(3) :: wj = (/delta,delta,gamma/)
25    CALL set_coefs_hevi(dt, &
26         (/ zero, (/2.*gamma,0.,0./), (/1-alpha,alpha,0./), wj /), &
27         (/ zero, (/gamma,gamma,0./), wj, wj /) )
28  END SUBROUTINE set_coefs_ark23
29
30  SUBROUTINE set_coefs_ark33(dt)
31    ! Fully-explicit RK3 scheme disguised as ARK
32    REAL(rstd) :: dt
[366]33    CALL set_coefs_rk(dt, (/ zero, (/.5,0.,0./), (/-1.,2.,0./), (/1./6.,2./3.,1./6./) /) )
[360]34  END SUBROUTINE set_coefs_ark33
35   
[366]36  SUBROUTINE set_coefs_rk(dt, ajl)
37    REAL(rstd) :: dt, ajl(3,4)
38    CALL set_coefs_hevi(dt,ajl,ajl)
39  END SUBROUTINE set_coefs_rk
40
[360]41  SUBROUTINE set_coefs_hevi(dt, ajl_slow, ajl_fast)
42    REAL(rstd) :: dt, ajl_slow(3,4), ajl_fast(3,4) ! fast/slow Butcher tableaus
43    INTEGER :: j
44    DO j=1,3
45       bjl(:,j) = dt*(ajl_slow(:,j+1)-ajl_slow(:,j))
46       cjl(:,j) = dt*(ajl_fast(:,j+1)-ajl_fast(:,j))
47       taujj(j) = dt*ajl_fast(j,j)
48    END DO
49    wj=dt*ajl_slow(:,4)
50  END SUBROUTINE set_coefs_hevi
51
52  SUBROUTINE HEVI_scheme(it, fluxt_zero)
53    USE time_mod
54    USE disvert_mod
[361]55    USE caldyn_hevi_mod
[480]56    USE omp_para
57    USE checksum_mod
[360]58    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
59    INTEGER :: it,j,l, ind
60    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
61
[529]62    CALL legacy_to_DEC(f_ps, f_u)
[360]63    DO j=1,nb_stage
[361]64       CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), &
[360]65            f_phis, f_ps,f_mass,f_theta_rhodz,f_u,f_q, &
[366]66            f_W, f_geopot, f_hflux, f_wflux, &
[360]67            f_dps_slow(:,j), f_dmass_slow(:,j), f_dtheta_rhodz_slow(:,j), &
[366]68            f_du_slow(:,j), f_du_fast(:,j), &
69            f_dPhi_slow(:,j), f_dPhi_fast(:,j), &
70            f_dW_slow(:,j), f_dW_fast(:,j) )
[361]71       ! accumulate mass fluxes for transport scheme
[480]72       
[360]73       DO ind=1,ndomain
74          IF (.NOT. assigned_domain(ind)) CYCLE
75          CALL swap_dimensions(ind)
76          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
77          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
78          CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, wj(j), fluxt_zero(ind))
79       END DO
80       ! update model state
81       DO l=1,j
82          IF(caldyn_eta==eta_mass) THEN
83             CALL update_2D(bjl(l,j), f_ps, f_dps_slow(:,l))
84          ELSE
[387]85             CALL update_3D(bjl(l,j), f_mass, f_dmass_slow(:,l))
[360]86          END IF
[387]87          CALL update_4D(bjl(l,j), f_theta_rhodz, f_dtheta_rhodz_slow(:,l))
88          CALL update_3D(bjl(l,j), f_u, f_du_slow(:,l))
89          CALL update_3D(cjl(l,j), f_u, f_du_fast(:,l))
[366]90          IF(.NOT. hydrostatic) THEN
[387]91             CALL update_3D(bjl(l,j), f_W, f_dW_slow(:,l))
92             CALL update_3D(cjl(l,j), f_W, f_dW_fast(:,l))
93             CALL update_3D(bjl(l,j), f_geopot, f_dPhi_slow(:,l))
94             CALL update_3D(cjl(l,j), f_geopot, f_dPhi_fast(:,l))
[366]95          END IF
[480]96!$OMP BARRIER       
[360]97       END DO
98    END DO
[529]99    CALL DEC_to_legacy(f_ps, f_u)
[360]100  END SUBROUTINE HEVI_scheme
101 
[387]102  SUBROUTINE update_4D(w, f_y, f_dy)
[360]103    USE dimensions
[387]104    USE grid_param, ONLY : nqdyn
[360]105    REAL(rstd) :: w
106    TYPE(t_field) :: f_y(:), f_dy(:)
[387]107    REAL(rstd), POINTER :: y(:,:,:), dy(:,:,:)
108    INTEGER :: ind, iq
109    IF(w /= 0.) THEN
110       DO ind=1,ndomain
111          IF (.NOT. assigned_domain(ind)) CYCLE
112          CALL swap_dimensions(ind)
113          dy=f_dy(ind); y=f_y(ind)
114          DO iq=1,nqdyn
115             CALL compute_update_3D(w,y(:,:,iq),dy(:,:,iq))
116          END DO
117       ENDDO
118    END IF
119  END SUBROUTINE update_4D
120   
121  SUBROUTINE update_3D(w, f_y, f_dy)
122    USE dimensions
123    REAL(rstd) :: w
124    TYPE(t_field) :: f_y(:), f_dy(:)
[360]125    REAL(rstd), POINTER :: y(:,:), dy(:,:)
126    INTEGER :: ind
[362]127    IF(w /= 0.) THEN
128       DO ind=1,ndomain
129          IF (.NOT. assigned_domain(ind)) CYCLE
130          CALL swap_dimensions(ind)
131          dy=f_dy(ind); y=f_y(ind)
[387]132          CALL compute_update_3D(w,y,dy)
[362]133       ENDDO
134    END IF
[387]135  END SUBROUTINE update_3D
[360]136   
[387]137  SUBROUTINE compute_update_3D(w, y, dy)
[360]138    USE omp_para
139    USE disvert_mod
140    REAL(rstd) :: w
141    REAL(rstd) :: y(:,:), dy(:,:)
142    INTENT(INOUT) :: y
143    INTENT(IN) :: dy
144    INTEGER :: l
145    DO l=ll_begin,ll_end
146       y(:,l)=y(:,l)+w*dy(:,l)
147    ENDDO
[387]148  END SUBROUTINE compute_update_3D
[360]149 
150  SUBROUTINE update_2D(w, f_y, f_dy)
[480]151  USE omp_para
[360]152    REAL(rstd) :: w
153    TYPE(t_field) :: f_y(:), f_dy(:)
154    REAL(rstd), POINTER :: y(:), dy(:)
155    INTEGER :: ind
156    DO ind=1,ndomain
157       IF (.NOT. assigned_domain(ind)) CYCLE
158       dy=f_dy(ind); y=f_y(ind)
[480]159       IF (is_omp_level_master) CALL compute_update_2D(w,y,dy)
[360]160    ENDDO
161  END SUBROUTINE update_2D
162   
163  SUBROUTINE compute_update_2D(w, y, dy)
164    REAL(rstd) :: w, y(:), dy(:)
165    INTENT(INOUT) :: y
166    INTENT(IN) :: dy
167    y(:)=y(:)+w*dy(:)
168  END SUBROUTINE compute_update_2D
169 
170END MODULE hevi_scheme_mod
Note: See TracBrowser for help on using the repository browser.