Context Navigation

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

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

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

Download in other formats: