Context Navigation

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

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_advect_tracer
	156	CALL init_check_conserve
[186]	157
[366]	158	CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
[818]	159
[820]	160	CALL init_guided(f_u,f_theta_rhodz)
[151]	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)
[162]	166	CALL init_message(f_mass,req_i0,req_mass0)
[151]	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)
[377]	170	CALL init_message(f_geopot,req_i0,req_geopot0)
	171	CALL init_message(f_W,req_i0,req_W0)
[151]	172
	173	END SUBROUTINE init_timeloop
[12]	174
[151]	175	SUBROUTINE timeloop
[360]	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
[584]	183	USE diagflux_mod
[360]	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
[528]	191	USE restart_mod
[360]	192	USE checksum_mod
	193	USE explicit_scheme_mod
	194	USE hevi_scheme_mod
	195	REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:)
[12]	196
[592]	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
[360]	200	LOGICAL, PARAMETER :: check_rhodz=.FALSE.
	201	INTEGER :: start_clock, stop_clock, rate_clock
[327]	202	LOGICAL,SAVE :: first_physic=.TRUE.
[360]	203	!$OMP THREADPRIVATE(first_physic)
	204
[295]	205	CALL switch_omp_distrib_level
[350]	206	CALL caldyn_BC(f_phis, f_geopot, f_wflux) ! set constant values in first/last interfaces
[132]	207
[360]	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.
[186]	224
[387]	225	IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
[592]	226	IF(diagflux_on) THEN
	227	adv_over_out = itau_adv*(1./itau_out)
	228	ELSE
	229	adv_over_out = 0.
	230	END IF
[377]	231
[360]	232	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)
[326]	233
[584]	234	Call trace_on
[186]	235
[413]	236	IF (xios_output) THEN ! we must call update_calendar before any XIOS output
[473]	237	IF (is_omp_master) CALL xios_update_calendar(1)
[413]	238	END IF
[528]	239
[569]	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)
[528]	242
[413]	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
[529]	245	!$OMP MASTER
[714]	246	CALL SYSTEM_CLOCK(start_clock, rate_clock)
[529]	247	!$OMP END MASTER
	248
[360]	249	DO it=itau0+1,itau0+itaumax
[364]	250	IF (is_master) CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock)
[413]	251
[714]	252	CALL enter_profile(id_timeloop)
	253
[364]	254	IF (xios_output) THEN
[473]	255
	256	IF(it>itau0+1 .AND. is_omp_master) CALL xios_update_calendar(it-itau0)
[364]	257	ELSE
	258	CALL update_time_counter(dt*it)
	259	END IF
[295]	260
[360]	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)
[377]	271	CALL wait_message(req_q0)
	272	IF(.NOT. hydrostatic) THEN
[499]	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)
[377]	277	END IF
[167]	278	ENDIF
[327]	279
[360]	280	CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
[326]	281
[714]	282	CALL enter_profile(id_dyn)
[360]	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
[714]	289	CALL exit_profile(id_dyn)
[594]	290
[714]	291	CALL enter_profile(id_dissip)
[360]	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)
[162]	303	END DO
[360]	304	ENDIF
	305
[714]	306	CALL enter_profile(id_diags)
[365]	307	CALL check_conserve_detailed(it, AAM_dyn, &
	308	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[714]	309	CALL exit_profile(id_diags)
[365]	310
[529]	311	CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du)
[360]	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
[162]	317	END IF
[365]	318
[714]	319	CALL enter_profile(id_diags)
[365]	320	CALL check_conserve_detailed(it, AAM_dissip, &
	321	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[714]	322	CALL exit_profile(id_diags)
[360]	323	END IF
[714]	324	CALL exit_profile(id_dissip)
[148]	325
[714]	326	CALL enter_profile(id_adv)
[360]	327	IF(MOD(it,itau_adv)==0) THEN
[592]	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
[594]	330	fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme
[595]	331	! At this point advect_tracer has obtained the halos of u and rhodz,
	332	! needed for correct computation of kinetic energy
[601]	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)
[595]	334
[594]	335	IF (check_rhodz) THEN ! check that rhodz is consistent with ps
[360]	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
[387]	344	IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
[159]	345	END IF
[714]	346	CALL exit_profile(id_adv)
[360]	347
[714]	348	CALL enter_profile(id_diags)
[360]	349	IF (MOD(it,itau_physics)==0) THEN
[365]	350	CALL check_conserve_detailed(it, AAM_dyn, &
	351	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[714]	352	CALL enter_profile(id_phys)
[741]	353	CALL physics(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
[714]	354	CALL exit_profile(id_phys)
[365]	355	CALL check_conserve_detailed(it, AAM_phys, &
	356	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[360]	357	!$OMP MASTER
	358	IF (first_physic) CALL SYSTEM_CLOCK(start_clock)
	359	!$OMP END MASTER
	360	first_physic=.FALSE.
[159]	361	END IF
[365]	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)
[413]	367	ENDIF
[407]	368
	369	IF (mod(it,itau_out)==0 ) THEN
	370	CALL transfert_request(f_u,req_e1_vect)
[413]	371	CALL write_output_fields_basic(.FALSE.,f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
[407]	372	ENDIF
[714]	373	CALL exit_profile(id_diags)
[407]	374
[714]	375	CALL exit_profile(id_timeloop)
[360]	376	END DO
	377
[569]	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)
[360]	380
[365]	381	CALL check_conserve_detailed(it, AAM_dyn, &
	382	f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
[360]	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
[364]	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
[365]	399	WRITE(,'(A,I7,A,F14.1)') "It No :",it," t :",dtit
[364]	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
[609]	405	total = per_step*itaumax
[364]	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
[714]	412	IF(MOD(it,itau_prof)==0) CALL print_profile
	413
[364]	414	END SUBROUTINE print_iteration
	415
[387]	416	SUBROUTINE copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
[377]	417	TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
[387]	418	REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
	419	INTEGER :: ind, iq
[377]	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)
[387]	427	DO iq=1, nqdyn
	428	q(:,:,iq) = theta_rhodz(:,:,iq)/rhodz(:,:)
	429	END DO
[377]	430	END DO
[387]	431	END SUBROUTINE copy_theta_to_q
[377]	432
[387]	433	SUBROUTINE copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
[377]	434	TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
[387]	435	REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
	436	INTEGER :: ind, iq
[377]	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)
[387]	444	DO iq=1,nqdyn
	445	theta_rhodz(:,:,iq) = rhodz(:,:)*q(:,:,iq)
	446	END DO
[377]	447	END DO
[387]	448	END SUBROUTINE copy_q_to_theta
[377]	449
[12]	450	END MODULE timeloop_gcm_mod

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

Download in other formats: