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

Last change on this file since 879 was 820, checked in by jisesh, 5 years ago

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