Context Navigation

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

Last change on this file since 732 was 732, checked in by dubos, 6 years ago
devel : more cleanup and reorganization in dynamics/
File size: 15.7 KB

Rev	Line
[12]	1	MODULE timeloop_gcm_mod
[714]	2	USE profiling_mod
[151]	3	USE icosa
[360]	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
[732]	10	USE caldyn_vars_mod
[350]	11	IMPLICIT NONE
[133]	12	PRIVATE
[12]	13
[186]	14	INTEGER, PARAMETER :: itau_sync=10
[377]	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
[714]	17	INTEGER :: itau_prof, id_timeloop, id_dyn, id_phys, id_dissip, id_adv, id_diags
[360]	18	PUBLIC :: init_timeloop, timeloop
[151]	19
[12]	20	CONTAINS
	21
[151]	22	SUBROUTINE init_timeloop
[360]	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
[12]	33
[360]	34	CHARACTER(len=255) :: def
[413]	35
[714]	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
[413]	43	CALL init_caldyn
[360]	44
[429]	45	! IF (xios_output) itau_out=1
[360]	46	IF (.NOT. enable_io) itau_out=HUGE(itau_out)
[377]	47
[714]	48	itau_prof=100
	49	CALL getin('itau_profiling',itau_prof)
	50
[377]	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
[376]	58	def='ARK2.3'
[360]	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
[159]	98	CALL allocate_field(f_phis,field_t,type_real,name='phis')
[360]	99	! Model state at current time step
[159]	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')
[360]	102	CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz')
[387]	103	CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,nqdyn,name='theta_rhodz')
[159]	104	CALL allocate_field(f_u,field_u,type_real,llm,name='u')
[366]	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')
[266]	107	CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q')
[360]	108	! Mass fluxes
[159]	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
[360]	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
[159]	117	CALL allocate_field(f_dps,field_t,type_real,name='dps')
[360]	118	CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
[387]	119	CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,nqdyn,name='dtheta_rhodz')
[360]	120	CALL allocate_field(f_du,field_u,type_real,llm,name='du')
	121	! Model state at previous time step (RK/MLF)
[159]	122	CALL allocate_field(f_psm1,field_t,type_real,name='psm1')
[162]	123	CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
[387]	124	CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,nqdyn,name='theta_rhodzm1')
[360]	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')
[387]	130	CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,nqdyn,name='dtheta_rhodz_fast')
[360]	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')
[366]	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')
[360]	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
[151]	141
[360]	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)
[387]	147	CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm,nqdyn)
[360]	148	CALL allocate_field(f_um2,field_u,type_real,llm)
[151]	149	END SELECT
	150
[295]	151	CALL init_theta2theta_rhodz
[151]	152	CALL init_dissip
[327]	153	CALL init_sponge
[354]	154	CALL init_observable
[151]	155	CALL init_guided
	156	CALL init_advect_tracer
	157	CALL init_check_conserve
[186]	158
[366]	159	CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
[151]	160
	161	CALL transfert_request(f_phis,req_i0)
	162	CALL transfert_request(f_phis,req_i1)
	163
	164	CALL init_message(f_ps,req_i0,req_ps0)
[162]	165	CALL init_message(f_mass,req_i0,req_mass0)
[151]	166	CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
	167	CALL init_message(f_u,req_e0_vect,req_u0)
	168	CALL init_message(f_q,req_i0,req_q0)
[377]	169	CALL init_message(f_geopot,req_i0,req_geopot0)
	170	CALL init_message(f_W,req_i0,req_W0)
[151]	171
	172	END SUBROUTINE init_timeloop
[12]	173
[151]	174	SUBROUTINE timeloop
[360]	175	USE dissip_gcm_mod
	176	USE sponge_mod
	177	USE observable_mod
	178	USE etat0_mod
	179	USE guided_mod
	180	USE caldyn_mod
	181	USE advect_tracer_mod
[584]	182	USE diagflux_mod
[360]	183	USE physics_mod
	184	USE mpipara
	185	USE transfert_mod
	186	USE check_conserve_mod
	187	USE xios_mod
	188	USE output_field_mod
	189	USE write_etat0_mod
[528]	190	USE restart_mod
[360]	191	USE checksum_mod
	192	USE explicit_scheme_mod
	193	USE hevi_scheme_mod
	194	REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:)
[12]	195
[592]	196	REAL(rstd) :: adv_over_out ! ratio itau_adv/itau_out
	197	INTEGER :: ind, it,i,j,l,n, stage
	198	LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating mass fluxes in time
[360]	199	LOGICAL, PARAMETER :: check_rhodz=.FALSE.
	200	INTEGER :: start_clock, stop_clock, rate_clock
[327]	201	LOGICAL,SAVE :: first_physic=.TRUE.
[360]	202	!$OMP THREADPRIVATE(first_physic)
	203
[295]	204	CALL switch_omp_distrib_level
[350]	205	CALL caldyn_BC(f_phis, f_geopot, f_wflux) ! set constant values in first/last interfaces
[132]	206
[360]	207	!$OMP BARRIER
	208	DO ind=1,ndomain
	209	IF (.NOT. assigned_domain(ind)) CYCLE
	210	CALL swap_dimensions(ind)
	211	CALL swap_geometry(ind)
	212	rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
	213	IF(caldyn_eta==eta_mass) THEN
	214	CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
	215	ELSE
	216	DO l=ll_begin,ll_end
	217	rhodz(:,l)=mass(:,l)
	218	ENDDO
	219	END IF
	220	END DO
	221	!$OMP BARRIER
	222	fluxt_zero=.TRUE.
[186]	223
[387]	224	IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
[592]	225	IF(diagflux_on) THEN
	226	adv_over_out = itau_adv*(1./itau_out)
	227	ELSE
	228	adv_over_out = 0.
	229	END IF
[377]	230
[360]	231	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)
[326]	232
[584]	233	Call trace_on
[186]	234
[413]	235	IF (xios_output) THEN ! we must call update_calendar before any XIOS output
[473]	236	IF (is_omp_master) CALL xios_update_calendar(1)
[413]	237	END IF
[528]	238
[569]	239	! IF (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)
	240	IF (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W)
[528]	241
[413]	242	CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
	243
[529]	244	!$OMP MASTER
[714]	245	CALL SYSTEM_CLOCK(start_clock, rate_clock)
[529]	246	!$OMP END MASTER
	247
[360]	248	DO it=itau0+1,itau0+itaumax
[364]	249	IF (is_master) CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock)
[413]	250
[714]	251	CALL enter_profile(id_timeloop)
	252
[364]	253	IF (xios_output) THEN
[473]	254
	255	IF(it>itau0+1 .AND. is_omp_master) CALL xios_update_calendar(it-itau0)
[364]	256	ELSE
	257	CALL update_time_counter(dt*it)
	258	END IF
[295]	259
[360]	260	IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN
	261	CALL send_message(f_ps,req_ps0)
	262	CALL wait_message(req_ps0)
	263	CALL send_message(f_mass,req_mass0)
	264	CALL wait_message(req_mass0)
	265	CALL send_message(f_theta_rhodz,req_theta_rhodz0)
	266	CALL wait_message(req_theta_rhodz0)
	267	CALL send_message(f_u,req_u0)
	268	CALL wait_message(req_u0)
	269	CALL send_message(f_q,req_q0)
[377]	270	CALL wait_message(req_q0)
	271	IF(.NOT. hydrostatic) THEN
[499]	272	CALL send_message(f_geopot,req_geopot0)
	273	CALL wait_message(req_geopot0)
	274	CALL send_message(f_W,req_W0)
	275	CALL wait_message(req_W0)
[377]	276	END IF
[167]	277	ENDIF
[327]	278
[360]	279	CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
[326]	280
[714]	281	CALL enter_profile(id_dyn)
[360]	282	SELECT CASE(scheme_family)
	283	CASE(explicit)
	284	CALL explicit_scheme(it, fluxt_zero)
	285	CASE(hevi)
	286	CALL HEVI_scheme(it, fluxt_zero)
	287	END SELECT
[714]	288	CALL exit_profile(id_dyn)
[594]	289
[714]	290	CALL enter_profile(id_dissip)
[360]	291	IF (MOD(it,itau_dissip)==0) THEN
	292
	293	IF(caldyn_eta==eta_mass) THEN
	294	!ym flush ps
	295	!$OMP BARRIER
	296	DO ind=1,ndomain
	297	IF (.NOT. assigned_domain(ind)) CYCLE
	298	CALL swap_dimensions(ind)
	299	CALL swap_geometry(ind)
	300	mass=f_mass(ind); ps=f_ps(ind);
	301	CALL compute_rhodz(.TRUE., ps, mass)
[162]	302	END DO
[360]	303	ENDIF
	304
[714]	305	CALL enter_profile(id_diags)
[365]	306	CALL check_conserve_detailed(it, AAM_dyn, &
	307	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[714]	308	CALL exit_profile(id_diags)
[365]	309
[529]	310	CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du)
[360]	311
	312	CALL euler_scheme(.FALSE.) ! update only u, theta
	313	IF (iflag_sponge > 0) THEN
	314	CALL sponge(f_u,f_du,f_theta_rhodz,f_dtheta_rhodz)
	315	CALL euler_scheme(.FALSE.) ! update only u, theta
[162]	316	END IF
[365]	317
[714]	318	CALL enter_profile(id_diags)
[365]	319	CALL check_conserve_detailed(it, AAM_dissip, &
	320	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[714]	321	CALL exit_profile(id_diags)
[360]	322	END IF
[714]	323	CALL exit_profile(id_dissip)
[148]	324
[714]	325	CALL enter_profile(id_adv)
[360]	326	IF(MOD(it,itau_adv)==0) THEN
[592]	327	CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step
	328	adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt) ! accumulate mass and fluxes if diagflux_on
[594]	329	fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme
[595]	330	! At this point advect_tracer has obtained the halos of u and rhodz,
	331	! needed for correct computation of kinetic energy
[601]	332	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)
[595]	333
[594]	334	IF (check_rhodz) THEN ! check that rhodz is consistent with ps
[360]	335	DO ind=1,ndomain
	336	IF (.NOT. assigned_domain(ind)) CYCLE
	337	CALL swap_dimensions(ind)
	338	CALL swap_geometry(ind)
	339	rhodz=f_rhodz(ind); ps=f_ps(ind);
	340	CALL compute_rhodz(.FALSE., ps, rhodz)
	341	END DO
	342	ENDIF
[387]	343	IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
[159]	344	END IF
[714]	345	CALL exit_profile(id_adv)
[360]	346
[714]	347	CALL enter_profile(id_diags)
[360]	348	IF (MOD(it,itau_physics)==0) THEN
[365]	349	CALL check_conserve_detailed(it, AAM_dyn, &
	350	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[714]	351	CALL enter_profile(id_phys)
[360]	352	CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
[714]	353	CALL exit_profile(id_phys)
[365]	354	CALL check_conserve_detailed(it, AAM_phys, &
	355	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[360]	356	!$OMP MASTER
	357	IF (first_physic) CALL SYSTEM_CLOCK(start_clock)
	358	!$OMP END MASTER
	359	first_physic=.FALSE.
[159]	360	END IF
[365]	361
	362	IF (MOD(it,itau_check_conserv)==0) THEN
	363	CALL check_conserve_detailed(it, AAM_dyn, &
	364	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
	365	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
[413]	366	ENDIF
[407]	367
	368	IF (mod(it,itau_out)==0 ) THEN
	369	CALL transfert_request(f_u,req_e1_vect)
[413]	370	CALL write_output_fields_basic(.FALSE.,f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
[407]	371	ENDIF
[714]	372	CALL exit_profile(id_diags)
[407]	373
[714]	374	CALL exit_profile(id_timeloop)
[360]	375	END DO
	376
[569]	377	! CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q)
	378	CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W)
[360]	379
[365]	380	CALL check_conserve_detailed(it, AAM_dyn, &
	381	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[360]	382	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
	383
	384	!$OMP MASTER
	385	CALL SYSTEM_CLOCK(stop_clock)
	386
	387	IF (mpi_rank==0) THEN
	388	PRINT ,"Time elapsed : ",(stop_clock-start_clock)1./rate_clock
	389	ENDIF
	390	!$OMP END MASTER
	391
	392	! CONTAINS
	393	END SUBROUTINE timeloop
[364]	394
	395	SUBROUTINE print_iteration(it,itau0,itaumax,start_clock,rate_clock)
	396	INTEGER :: it, itau0, itaumax, start_clock, stop_clock, rate_clock, throughput
	397	REAL :: per_step,total, elapsed
[365]	398	WRITE(,'(A,I7,A,F14.1)') "It No :",it," t :",dtit
[364]	399	IF(MOD(it,10)==0) THEN
	400	CALL SYSTEM_CLOCK(stop_clock)
	401	elapsed = (stop_clock-start_clock)*1./rate_clock
	402	per_step = elapsed/(it-itau0)
	403	throughput = dt/per_step
[609]	404	total = per_step*itaumax
[364]	405	WRITE(*,'(A,I5,A,F6.2,A,I6)') 'Time spent (s):',INT(elapsed), &
	406	' -- ms/step : ', 1000*per_step, &
	407	' -- Throughput :', throughput
	408	WRITE(*,'(A,I5,A,I5)') 'Whole job (min) :', INT(total/60.), &
	409	' -- Completion in (min) : ', INT((total-elapsed)/60.)
	410	END IF
[714]	411	IF(MOD(it,itau_prof)==0) CALL print_profile
	412
[364]	413	END SUBROUTINE print_iteration
	414
[387]	415	SUBROUTINE copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
[377]	416	TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
[387]	417	REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
	418	INTEGER :: ind, iq
[377]	419	DO ind=1, ndomain
	420	IF (.NOT. assigned_domain(ind)) CYCLE
	421	CALL swap_dimensions(ind)
	422	CALL swap_geometry(ind)
	423	theta_rhodz=f_theta_rhodz(ind)
	424	rhodz=f_rhodz(ind)
	425	q=f_q(ind)
[387]	426	DO iq=1, nqdyn
	427	q(:,:,iq) = theta_rhodz(:,:,iq)/rhodz(:,:)
	428	END DO
[377]	429	END DO
[387]	430	END SUBROUTINE copy_theta_to_q
[377]	431
[387]	432	SUBROUTINE copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
[377]	433	TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
[387]	434	REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
	435	INTEGER :: ind, iq
[377]	436	DO ind=1, ndomain
	437	IF (.NOT. assigned_domain(ind)) CYCLE
	438	CALL swap_dimensions(ind)
	439	CALL swap_geometry(ind)
	440	theta_rhodz=f_theta_rhodz(ind)
	441	rhodz=f_rhodz(ind)
	442	q=f_q(ind)
[387]	443	DO iq=1,nqdyn
	444	theta_rhodz(:,:,iq) = rhodz(:,:)*q(:,:,iq)
	445	END DO
[377]	446	END DO
[387]	447	END SUBROUTINE copy_q_to_theta
[377]	448
[12]	449	END MODULE timeloop_gcm_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats: