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

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

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