source: codes/icosagcm/trunk/src/timeloop_gcm.f90 @ 376

Last change on this file since 376 was 376, checked in by dubos, 8 years ago

Bugfix : define again ARK2.3 as default time scheme
File size: 11.7 KB
Line 
1MODULE timeloop_gcm_mod
2  USE icosa
3  USE disvert_mod
4  USE trace
5  USE omp_para
6  USE euler_scheme_mod
7  USE explicit_scheme_mod
8  USE hevi_scheme_mod
9  IMPLICIT NONE
10  PRIVATE
11 
12  INTEGER, PARAMETER :: itau_sync=10
13  TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0
14
15  PUBLIC :: init_timeloop, timeloop
16
17CONTAINS
18 
19  SUBROUTINE init_timeloop
20    USE dissip_gcm_mod
21    USE observable_mod
22    USE caldyn_mod
23    USE etat0_mod
24    USE guided_mod
25    USE advect_tracer_mod
26    USE check_conserve_mod
27    USE output_field_mod
28    USE theta2theta_rhodz_mod
29    USE sponge_mod
30
31    CHARACTER(len=255) :: def
32   
33    IF (xios_output) itau_out=1
34    IF (.NOT. enable_io) itau_out=HUGE(itau_out)
35   
36    def='ARK2.3'
37    CALL getin('time_scheme',def)
38    SELECT CASE (TRIM(def))
39    CASE('euler')
40       scheme_family=explicit
41       scheme=euler
42       nb_stage=1
43    CASE ('runge_kutta')
44       scheme_family=explicit
45       scheme=rk4
46       nb_stage=4
47    CASE ('RK2.5')
48       scheme_family=explicit
49       scheme=rk25
50       nb_stage=5
51    CASE ('leapfrog_matsuno')
52       scheme_family=explicit
53       scheme=mlf
54       matsuno_period=5
55       CALL getin('matsuno_period',matsuno_period)
56       nb_stage=matsuno_period+1
57    CASE('ARK2.3')
58       scheme_family=hevi
59       scheme=ark23
60       nb_stage=3
61       CALL set_coefs_ark23(dt)
62    CASE('ARK3.3')
63       scheme_family=hevi
64       scheme=ark33
65       nb_stage=3
66       CALL set_coefs_ark33(dt)
67    CASE ('none')
68       nb_stage=0
69    CASE default
70       PRINT*,'Bad selector for variable scheme : <', TRIM(def),             &
71            ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>,<RK2.5>,<ARK2.3>'
72       STOP
73    END SELECT
74   
75    ! Time-independant orography
76    CALL allocate_field(f_phis,field_t,type_real,name='phis')
77    ! Model state at current time step
78    CALL allocate_field(f_ps,field_t,type_real, name='ps')
79    CALL allocate_field(f_mass,field_t,type_real,llm,name='mass')
80    CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz')
81    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz')
82    CALL allocate_field(f_u,field_u,type_real,llm,name='u')
83    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')
84    CALL allocate_field(f_W,field_t,type_real,llm+1,name='W')
85    CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q')
86    ! Mass fluxes
87    CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
88    CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
89    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes
90    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt')
91   
92    SELECT CASE(scheme_family)
93    CASE(explicit)
94       ! Trends
95       CALL allocate_field(f_dps,field_t,type_real,name='dps')
96       CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
97       CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz')
98       CALL allocate_field(f_du,field_u,type_real,llm,name='du')
99       ! Model state at previous time step (RK/MLF)
100       CALL allocate_field(f_psm1,field_t,type_real,name='psm1')
101       CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
102       CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1')
103       CALL allocate_field(f_um1,field_u,type_real,llm,name='um1')
104    CASE(hevi)
105       ! Trends
106       CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow')
107       CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow')
108       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,name='dtheta_rhodz_fast')
109       CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow')
110       CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast')
111       CALL allocate_fields(nb_stage,f_dW_slow, field_t,type_real,llm+1,name='dW_slow')
112       CALL allocate_fields(nb_stage,f_dW_fast, field_t,type_real,llm+1,name='dW_fast')
113       CALL allocate_fields(nb_stage,f_dPhi_slow, field_t,type_real,llm+1,name='dPhi_slow')
114       CALL allocate_fields(nb_stage,f_dPhi_fast, field_t,type_real,llm+1,name='dPhi_fast')
115       f_dps => f_dps_slow(:,1)
116       f_du => f_du_slow(:,1)
117       f_dtheta_rhodz => f_dtheta_rhodz_slow(:,1)
118    END SELECT
119
120    SELECT CASE(scheme)
121    CASE(mlf)
122       ! Model state 2 time steps ago (MLF)
123       CALL allocate_field(f_psm2,field_t,type_real)
124       CALL allocate_field(f_massm2,field_t,type_real,llm)
125       CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
126       CALL allocate_field(f_um2,field_u,type_real,llm)
127    END SELECT
128
129    CALL init_theta2theta_rhodz
130    CALL init_dissip
131    CALL init_sponge
132    CALL init_observable
133    CALL init_caldyn
134    CALL init_guided
135    CALL init_advect_tracer
136    CALL init_check_conserve
137
138    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
139
140    CALL transfert_request(f_phis,req_i0) 
141    CALL transfert_request(f_phis,req_i1) 
142    CALL writefield("phis",f_phis,once=.TRUE.)
143
144    CALL init_message(f_ps,req_i0,req_ps0)
145    CALL init_message(f_mass,req_i0,req_mass0)
146    CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
147    CALL init_message(f_u,req_e0_vect,req_u0)
148    CALL init_message(f_q,req_i0,req_q0)
149
150  END SUBROUTINE init_timeloop
151 
152  SUBROUTINE timeloop
153    USE dissip_gcm_mod
154    USE sponge_mod
155    USE observable_mod
156    USE etat0_mod
157    USE guided_mod
158    USE caldyn_mod
159    USE advect_tracer_mod
160    USE physics_mod
161    USE mpipara
162    USE transfert_mod
163    USE check_conserve_mod
164    USE xios_mod
165    USE output_field_mod
166    USE write_etat0_mod
167    USE checksum_mod
168    USE explicit_scheme_mod
169    USE hevi_scheme_mod
170    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:)
171
172    INTEGER :: ind
173    INTEGER :: it,i,j,l,n,  stage
174    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
175    LOGICAL, PARAMETER :: check_rhodz=.FALSE.
176    INTEGER :: start_clock, stop_clock, rate_clock
177    LOGICAL,SAVE :: first_physic=.TRUE.
178    !$OMP THREADPRIVATE(first_physic)   
179
180    CALL switch_omp_distrib_level
181    CALL caldyn_BC(f_phis, f_geopot, f_wflux) ! set constant values in first/last interfaces
182
183    !$OMP BARRIER
184    DO ind=1,ndomain
185       IF (.NOT. assigned_domain(ind)) CYCLE
186       CALL swap_dimensions(ind)
187       CALL swap_geometry(ind)
188       rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
189       IF(caldyn_eta==eta_mass) THEN
190          CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
191       ELSE
192          DO l=ll_begin,ll_end
193             rhodz(:,l)=mass(:,l)
194          ENDDO
195       END IF
196    END DO
197    !$OMP BARRIER
198    fluxt_zero=.TRUE.
199
200    !$OMP MASTER
201    CALL SYSTEM_CLOCK(start_clock)
202    CALL SYSTEM_CLOCK(count_rate=rate_clock)
203    !$OMP END MASTER   
204
205    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0) 
206
207    CALL trace_on
208
209    DO it=itau0+1,itau0+itaumax
210
211       IF (is_master) CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock)
212       IF (xios_output) THEN
213          CALL xios_update_calendar(it)
214       ELSE
215          CALL update_time_counter(dt*it)
216       END IF
217
218       IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN
219          CALL send_message(f_ps,req_ps0)
220          CALL wait_message(req_ps0)
221          CALL send_message(f_mass,req_mass0)
222          CALL wait_message(req_mass0)
223          CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
224          CALL wait_message(req_theta_rhodz0) 
225          CALL send_message(f_u,req_u0)
226          CALL wait_message(req_u0)
227          CALL send_message(f_q,req_q0) 
228          CALL wait_message(req_q0) 
229       ENDIF
230
231       IF (mod(it,itau_out)==0 ) THEN
232          CALL transfert_request(f_u,req_e1_vect)
233          CALL write_output_fields_basic(f_ps, f_mass, f_geopot, f_u, f_W, f_q)
234       ENDIF
235
236       CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
237
238       SELECT CASE(scheme_family)
239       CASE(explicit)
240          CALL explicit_scheme(it, fluxt_zero)
241       CASE(hevi)
242          CALL HEVI_scheme(it, fluxt_zero)
243       END SELECT
244       
245       IF (MOD(it,itau_dissip)==0) THEN
246         
247          IF(caldyn_eta==eta_mass) THEN
248             !ym flush ps
249             !$OMP BARRIER
250             DO ind=1,ndomain
251                IF (.NOT. assigned_domain(ind)) CYCLE
252                CALL swap_dimensions(ind)
253                CALL swap_geometry(ind)
254                mass=f_mass(ind); ps=f_ps(ind);
255                CALL compute_rhodz(.TRUE., ps, mass)
256             END DO
257          ENDIF
258         
259          CALL check_conserve_detailed(it, AAM_dyn, &
260               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
261       
262          CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz)
263         
264          CALL euler_scheme(.FALSE.)  ! update only u, theta
265          IF (iflag_sponge > 0) THEN
266             CALL sponge(f_u,f_du,f_theta_rhodz,f_dtheta_rhodz)
267             CALL euler_scheme(.FALSE.)  ! update only u, theta
268          END IF
269
270          CALL check_conserve_detailed(it, AAM_dissip, &
271               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
272       END IF
273       
274       IF(MOD(it,itau_adv)==0) THEN
275          CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
276          fluxt_zero=.TRUE.
277          ! FIXME : check that rhodz is consistent with ps
278          IF (check_rhodz) THEN
279             DO ind=1,ndomain
280                IF (.NOT. assigned_domain(ind)) CYCLE
281                CALL swap_dimensions(ind)
282                CALL swap_geometry(ind)
283                rhodz=f_rhodz(ind); ps=f_ps(ind);
284                CALL compute_rhodz(.FALSE., ps, rhodz)   
285             END DO
286          ENDIF
287       END IF
288       
289       IF (MOD(it,itau_physics)==0) THEN
290          CALL check_conserve_detailed(it, AAM_dyn, &
291            f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
292          CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)       
293          CALL check_conserve_detailed(it, AAM_phys, &
294               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
295          !$OMP MASTER
296          IF (first_physic) CALL SYSTEM_CLOCK(start_clock)
297          !$OMP END MASTER   
298          first_physic=.FALSE.
299       END IF
300
301       IF (MOD(it,itau_check_conserv)==0) THEN
302          CALL check_conserve_detailed(it, AAM_dyn, &
303               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
304          CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
305       ENDIF       
306    END DO
307   
308    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q) 
309   
310    CALL check_conserve_detailed(it, AAM_dyn, &
311         f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
312    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
313   
314    !$OMP MASTER
315    CALL SYSTEM_CLOCK(stop_clock)
316   
317    IF (mpi_rank==0) THEN
318       PRINT *,"Time elapsed : ",(stop_clock-start_clock)*1./rate_clock 
319    ENDIF
320    !$OMP END MASTER
321   
322    ! CONTAINS
323  END SUBROUTINE timeloop
324
325  SUBROUTINE print_iteration(it,itau0,itaumax,start_clock,rate_clock)
326    INTEGER :: it, itau0, itaumax, start_clock, stop_clock, rate_clock, throughput
327    REAL :: per_step,total, elapsed
328    WRITE(*,'(A,I7,A,F14.1)') "It No :",it,"   t :",dt*it
329    IF(MOD(it,10)==0) THEN
330       CALL SYSTEM_CLOCK(stop_clock)
331       elapsed = (stop_clock-start_clock)*1./rate_clock
332       per_step = elapsed/(it-itau0)
333       throughput = dt/per_step
334       total = per_step*(itaumax-itau0)
335       WRITE(*,'(A,I5,A,F6.2,A,I6)') 'Time spent (s):',INT(elapsed), &
336            '  -- ms/step : ', 1000*per_step, &
337            '  -- Throughput :', throughput
338       WRITE(*,'(A,I5,A,I5)') 'Whole job (min) :', INT(total/60.), &
339            '  -- Completion in (min) : ', INT((total-elapsed)/60.)
340    END IF
341  END SUBROUTINE print_iteration
342
343END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.