source: codes/icosagcm/devel/src/time/timeloop_gcm.f90 @ 595

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

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