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

Last change on this file since 915 was 909, checked in by dubos, 5 years ago

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