Context Navigation

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

Last change on this file since 318 was 295, checked in by ymipsl, 10 years ago

Merging OpenMP parallisme mode : by subdomain and on vertical level.
This feature is actually experimental but may be retro-compatible with the last method based only on subdomain

YM

File size: 19.4 KB

Rev	Line
[12]	1	MODULE timeloop_gcm_mod
[151]	2	USE transfert_mod
	3	USE icosa
[133]	4	PRIVATE
[12]	5
[151]	6	PUBLIC :: init_timeloop, timeloop
[133]	7
	8	INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
[186]	9	INTEGER, PARAMETER :: itau_sync=10
[133]	10
[186]	11	TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0
[151]	12
[186]	13	TYPE(t_field),POINTER,SAVE :: f_q(:)
	14	TYPE(t_field),POINTER,SAVE :: f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:)
	15	TYPE(t_field),POINTER,SAVE :: f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:)
	16	TYPE(t_field),POINTER,SAVE :: f_u(:),f_um1(:),f_um2(:), f_du(:)
	17	TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:)
	18	TYPE(t_field),POINTER,SAVE :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)
[151]	19
[186]	20	INTEGER,SAVE :: nb_stage, matsuno_period, scheme
	21	!$OMP THREADPRIVATE(nb_stage, matsuno_period, scheme)
[151]	22
[12]	23	CONTAINS
	24
[151]	25	SUBROUTINE init_timeloop
[19]	26	USE icosa
[15]	27	USE dissip_gcm_mod
[17]	28	USE caldyn_mod
[12]	29	USE etat0_mod
[159]	30	USE disvert_mod
[17]	31	USE guided_mod
	32	USE advect_tracer_mod
[81]	33	USE physics_mod
[131]	34	USE mpipara
[151]	35	USE omp_para
[145]	36	USE trace
[148]	37	USE transfert_mod
[151]	38	USE check_conserve_mod
[171]	39	USE output_field_mod
	40	USE write_field
[295]	41	USE theta2theta_rhodz_mod
[12]	42	IMPLICIT NONE
	43
[159]	44	CHARACTER(len=255) :: def
[17]	45
[149]	46
[171]	47	IF (xios_output) itau_out=1
[186]	48	IF (.NOT. enable_io) itau_out=HUGE(itau_out)
[129]	49
[159]	50	! Time-independant orography
	51	CALL allocate_field(f_phis,field_t,type_real,name='phis')
	52	! Trends
	53	CALL allocate_field(f_du,field_u,type_real,llm,name='du')
	54	CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz')
	55	! Model state at current time step (RK/MLF/Euler)
	56	CALL allocate_field(f_ps,field_t,type_real, name='ps')
	57	CALL allocate_field(f_mass,field_t,type_real,llm,name='mass')
	58	CALL allocate_field(f_u,field_u,type_real,llm,name='u')
	59	CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz')
	60	! Model state at previous time step (RK/MLF)
	61	CALL allocate_field(f_um1,field_u,type_real,llm,name='um1')
	62	CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1')
	63	! Tracers
[266]	64	CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q')
[159]	65	CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz')
	66	! Mass fluxes
	67	CALL allocate_field(f_hflux,field_u,type_real,llm) ! instantaneous mass fluxes
	68	CALL allocate_field(f_hfluxt,field_u,type_real,llm) ! mass "fluxes" accumulated in time
	69	CALL allocate_field(f_wflux,field_t,type_real,llm+1) ! vertical mass fluxes
[162]	70	CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
[151]	71
[159]	72	IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default)
	73	CALL allocate_field(f_dps,field_t,type_real,name='dps')
	74	CALL allocate_field(f_psm1,field_t,type_real,name='psm1')
	75	CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt')
	76	! the following are unused but must point to something
[162]	77	! f_massm1 => f_mass
[159]	78	ELSE ! eta = Lagrangian vertical coordinate
[162]	79	CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
[159]	80	! the following are unused but must point to something
	81	f_wfluxt => f_wflux
	82	f_dps => f_phis
	83	f_psm1 => f_phis
	84	END IF
[151]	85
[159]	86	def='runge_kutta'
	87	CALL getin('scheme',def)
	88
	89	SELECT CASE (TRIM(def))
[151]	90	CASE('euler')
	91	scheme=euler
	92	nb_stage=1
	93	CASE ('runge_kutta')
	94	scheme=rk4
	95	nb_stage=4
	96	CASE ('leapfrog_matsuno')
	97	scheme=mlf
	98	matsuno_period=5
	99	CALL getin('matsuno_period',matsuno_period)
	100	nb_stage=matsuno_period+1
[129]	101	! Model state 2 time steps ago (MLF)
[151]	102	CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
	103	CALL allocate_field(f_um2,field_u,type_real,llm)
[159]	104	IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default)
	105	CALL allocate_field(f_psm2,field_t,type_real)
	106	! the following are unused but must point to something
	107	f_massm2 => f_mass
	108	ELSE ! eta = Lagrangian vertical coordinate
	109	CALL allocate_field(f_massm2,field_t,type_real,llm)
	110	! the following are unused but must point to something
	111	f_psm2 => f_phis
	112	END IF
	113
[151]	114	CASE default
[159]	115	PRINT*,'Bad selector for variable scheme : <', TRIM(def), &
[151]	116	' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
	117	STOP
	118	END SELECT
	119
[295]	120	CALL init_theta2theta_rhodz
[151]	121	CALL init_dissip
	122	CALL init_caldyn
	123	CALL init_guided
	124	CALL init_advect_tracer
	125	CALL init_check_conserve
	126	CALL init_physics
[295]	127
[186]	128
[159]	129	CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
[151]	130
	131	CALL transfert_request(f_phis,req_i0)
	132	CALL transfert_request(f_phis,req_i1)
	133	CALL writefield("phis",f_phis,once=.TRUE.)
	134
	135	CALL init_message(f_ps,req_i0,req_ps0)
[162]	136	CALL init_message(f_mass,req_i0,req_mass0)
[151]	137	CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
	138	CALL init_message(f_u,req_e0_vect,req_u0)
	139	CALL init_message(f_q,req_i0,req_q0)
	140
	141	END SUBROUTINE init_timeloop
[12]	142
[151]	143	SUBROUTINE timeloop
	144	USE icosa
	145	USE dissip_gcm_mod
[159]	146	USE disvert_mod
[151]	147	USE caldyn_mod
[162]	148	USE caldyn_gcm_mod, ONLY : req_ps, req_mass
[151]	149	USE etat0_mod
	150	USE guided_mod
	151	USE advect_tracer_mod
	152	USE physics_mod
	153	USE mpipara
	154	USE omp_para
	155	USE trace
	156	USE transfert_mod
	157	USE check_conserve_mod
[171]	158	USE xios_mod
	159	USE output_field_mod
[266]	160	USE write_etat0_mod
[295]	161	USE checksum_mod
[151]	162	IMPLICIT NONE
	163	REAL(rstd),POINTER :: q(:,:,:)
[159]	164	REAL(rstd),POINTER :: phis(:), ps(:) ,psm1(:), psm2(:), dps(:)
	165	REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:)
	166	REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:)
	167	REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:)
[151]	168	REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
[12]	169
[151]	170	INTEGER :: ind
	171	INTEGER :: it,i,j,n, stage
	172	LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
	173	LOGICAL, PARAMETER :: check=.FALSE.
[186]	174	INTEGER :: start_clock
	175	INTEGER :: stop_clock
	176	INTEGER :: rate_clock
[295]	177	INTEGER :: l
[186]	178
[266]	179
	180	! CALL write_etat0(f_ps, f_phis,f_theta_rhodz,f_u,f_q)
	181	! CALL read_start(f_ps,f_mass,f_phis,f_theta_rhodz,f_u,f_q)
	182	! CALL write_restart(f_ps,f_mass,f_phis,f_theta_rhodz,f_u,f_q)
	183
[295]	184	CALL switch_omp_distrib_level
[159]	185	CALL caldyn_BC(f_phis, f_wflux) ! set constant values in first/last interfaces
[295]	186
	187	!$OMP BARRIER
[133]	188	DO ind=1,ndomain
[202]	189	IF (.NOT. assigned_domain(ind)) CYCLE
[133]	190	CALL swap_dimensions(ind)
	191	CALL swap_geometry(ind)
[162]	192	rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
	193	IF(caldyn_eta==eta_mass) THEN
	194	CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
	195	ELSE
[295]	196	DO l=ll_begin,ll_end
	197	rhodz(:,l)=mass(:,l)
	198	ENDDO
[162]	199	END IF
[133]	200	END DO
[295]	201	!$OMP BARRIER
[138]	202	fluxt_zero=.TRUE.
[132]	203
[186]	204	!$OMP MASTER
	205	CALL SYSTEM_CLOCK(start_clock)
	206	!$OMP END MASTER
	207
[266]	208	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)
	209
[186]	210	CALL trace_on
	211
[266]	212	DO it=itau0+1,itau0+itaumax
[171]	213
[200]	214	IF (xios_output) CALL xios_update_calendar(it)
[266]	215	IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN
[151]	216	CALL send_message(f_ps,req_ps0)
[186]	217	CALL wait_message(req_ps0)
[162]	218	CALL send_message(f_mass,req_mass0)
[186]	219	CALL wait_message(req_mass0)
[151]	220	CALL send_message(f_theta_rhodz,req_theta_rhodz0)
[186]	221	CALL wait_message(req_theta_rhodz0)
[151]	222	CALL send_message(f_u,req_u0)
[186]	223	CALL wait_message(req_u0)
[151]	224	CALL send_message(f_q,req_q0)
	225	CALL wait_message(req_q0)
[186]	226
	227	! CALL wait_message(req_ps0)
	228	! CALL wait_message(req_mass0)
	229	! CALL wait_message(req_theta_rhodz0)
	230	! CALL wait_message(req_u0)
	231	! CALL wait_message(req_q0)
[148]	232	ENDIF
[186]	233
[295]	234	IF (is_master) PRINT ,"It No :",It," t :",dtIt
	235
[63]	236	IF (mod(it,itau_out)==0 ) THEN
[81]	237	CALL update_time_counter(dt*it)
[171]	238	CALL output_field("q",f_q)
[63]	239	ENDIF
[151]	240
	241	CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
[129]	242
	243	DO stage=1,nb_stage
	244	CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
[159]	245	f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, &
[162]	246	f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
[133]	247	SELECT CASE (scheme)
	248	CASE(euler)
	249	CALL euler_scheme(.TRUE.)
	250	CASE (rk4)
[129]	251	CALL rk_scheme(stage)
[133]	252	CASE (mlf)
[129]	253	CALL leapfrog_matsuno_scheme(stage)
[133]	254	CASE DEFAULT
[129]	255	STOP
	256	END SELECT
	257	END DO
[130]	258
[266]	259	IF (MOD(it,itau_dissip)==0) THEN
[186]	260	! CALL send_message(f_ps,req_ps)
	261	! CALL wait_message(req_ps)
	262
[167]	263	IF(caldyn_eta==eta_mass) THEN
[295]	264	!ym flush ps
	265	!$OMP BARRIER
[167]	266	DO ind=1,ndomain
[186]	267	IF (.NOT. assigned_domain(ind)) CYCLE
[167]	268	CALL swap_dimensions(ind)
	269	CALL swap_geometry(ind)
	270	mass=f_mass(ind); ps=f_ps(ind);
	271	CALL compute_rhodz(.TRUE., ps, mass)
	272	END DO
	273	ENDIF
[186]	274	! CALL send_message(f_mass,req_mass)
	275	! CALL wait_message(req_mass)
[167]	276	CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz)
[295]	277
[186]	278	! CALL send_message(f_mass,req_mass)
	279	! CALL wait_message(req_mass)
[167]	280	CALL euler_scheme(.FALSE.) ! update only u, theta
	281	END IF
[130]	282
[266]	283	IF(MOD(it,itau_adv)==0) THEN
[138]	284
[135]	285	CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step
[134]	286	fluxt_zero=.TRUE.
[138]	287
	288	! FIXME : check that rhodz is consistent with ps
[148]	289	IF (check) THEN
	290	DO ind=1,ndomain
[186]	291	IF (.NOT. assigned_domain(ind)) CYCLE
[148]	292	CALL swap_dimensions(ind)
	293	CALL swap_geometry(ind)
[151]	294	rhodz=f_rhodz(ind); ps=f_ps(ind);
[148]	295	CALL compute_rhodz(.FALSE., ps, rhodz)
	296	END DO
	297	ENDIF
[151]	298
[133]	299	END IF
[151]	300
[281]	301	CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
[295]	302
	303	IF (MOD(it,itau_check_conserv)==0) THEN
	304	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
	305	ENDIF
[266]	306
	307	ENDDO
[151]	308
[266]	309	CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q)
[186]	310
[266]	311	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
[151]	312
[186]	313	!$OMP MASTER
[266]	314	CALL SYSTEM_CLOCK(stop_clock)
	315	CALL SYSTEM_CLOCK(count_rate=rate_clock)
[186]	316
[266]	317	IF (mpi_rank==0) THEN
	318	PRINT ,"Time elapsed : ",(stop_clock-start_clock)1./rate_clock
	319	ENDIF
[186]	320	!$OMP END MASTER
[129]	321
[266]	322	CONTAINS
[12]	323
[130]	324	SUBROUTINE Euler_scheme(with_dps)
[12]	325	IMPLICIT NONE
[130]	326	LOGICAL :: with_dps
	327	INTEGER :: ind
[148]	328	INTEGER :: i,j,ij,l
[145]	329	CALL trace_start("Euler_scheme")
	330
[130]	331	DO ind=1,ndomain
[186]	332	IF (.NOT. assigned_domain(ind)) CYCLE
[138]	333	CALL swap_dimensions(ind)
	334	CALL swap_geometry(ind)
[148]	335
[162]	336	IF(with_dps) THEN ! update ps/mass
	337	IF(caldyn_eta==eta_mass) THEN ! update ps
	338	ps=f_ps(ind) ; dps=f_dps(ind) ;
[295]	339	IF (is_omp_first_level) THEN
[174]	340	!$SIMD
	341	DO ij=ij_begin,ij_end
	342	ps(ij)=ps(ij)+dt*dps(ij)
[162]	343	ENDDO
	344	ENDIF
	345	ELSE ! update mass
	346	mass=f_mass(ind) ; dmass=f_dmass(ind) ;
	347	DO l=1,llm
[174]	348	!$SIMD
	349	DO ij=ij_begin,ij_end
	350	mass(ij,l)=mass(ij,l)+dt*dmass(ij,l)
[162]	351	ENDDO
	352	END DO
	353	END IF
	354
	355	hflux=f_hflux(ind); hfluxt=f_hfluxt(ind)
	356	wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
	357	CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
	358	END IF ! update ps/mass
[148]	359
[130]	360	u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
	361	du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
[148]	362
[151]	363	DO l=ll_begin,ll_end
[174]	364	!$SIMD
	365	DO ij=ij_begin,ij_end
[148]	366	u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l)
	367	u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
	368	u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
	369	theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l)
	370	ENDDO
	371	ENDDO
[130]	372	ENDDO
[133]	373
[145]	374	CALL trace_end("Euler_scheme")
	375
[12]	376	END SUBROUTINE Euler_scheme
[120]	377
[129]	378	SUBROUTINE RK_scheme(stage)
[120]	379	IMPLICIT NONE
	380	INTEGER :: ind, stage
[129]	381	REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
[120]	382	REAL(rstd) :: tau
[148]	383	INTEGER :: i,j,ij,l
[145]	384
	385	CALL trace_start("RK_scheme")
[120]	386
	387	tau = dt*coef(stage)
[151]	388
[159]	389	! if mass coordinate, deal with ps first on one core
	390	IF(caldyn_eta==eta_mass) THEN
[295]	391	IF (is_omp_first_level) THEN
[162]	392
[159]	393	DO ind=1,ndomain
[186]	394	IF (.NOT. assigned_domain(ind)) CYCLE
[159]	395	CALL swap_dimensions(ind)
	396	CALL swap_geometry(ind)
[162]	397	ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind)
[159]	398
	399	IF (stage==1) THEN ! first stage : save model state in XXm1
[174]	400	!$SIMD
	401	DO ij=ij_begin,ij_end
	402	psm1(ij)=ps(ij)
	403	ENDDO
[159]	404	ENDIF
	405
	406	! updates are of the form x1 := x0 + tau*f(x1)
[174]	407	!$SIMD
	408	DO ij=ij_begin,ij_end
	409	ps(ij)=psm1(ij)+tau*dps(ij)
[151]	410	ENDDO
[159]	411	ENDDO
	412	ENDIF
[186]	413	! CALL send_message(f_ps,req_ps)
	414	!ym no overlap for now
	415	! CALL wait_message(req_ps)
[162]	416
	417	ELSE ! Lagrangian coordinate, deal with mass
	418	DO ind=1,ndomain
[186]	419	IF (.NOT. assigned_domain(ind)) CYCLE
[162]	420	CALL swap_dimensions(ind)
	421	CALL swap_geometry(ind)
	422	mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind)
	423
	424	IF (stage==1) THEN ! first stage : save model state in XXm1
	425	DO l=ll_begin,ll_end
[174]	426	!$SIMD
	427	DO ij=ij_begin,ij_end
	428	massm1(ij,l)=mass(ij,l)
	429	ENDDO
[162]	430	ENDDO
	431	END IF
	432
	433	! updates are of the form x1 := x0 + tau*f(x1)
	434	DO l=ll_begin,ll_end
[174]	435	!$SIMD
	436	DO ij=ij_begin,ij_end
	437	mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
[162]	438	ENDDO
	439	ENDDO
	440	END DO
[186]	441	! CALL send_message(f_mass,req_mass)
	442	!ym no overlap for now
	443	! CALL wait_message(req_mass)
[162]	444
[159]	445	END IF
[151]	446
[159]	447	! now deal with other prognostic variables
[120]	448	DO ind=1,ndomain
[186]	449	IF (.NOT. assigned_domain(ind)) CYCLE
[138]	450	CALL swap_dimensions(ind)
	451	CALL swap_geometry(ind)
[162]	452	u=f_u(ind) ; du=f_du(ind) ; um1=f_um1(ind)
	453	theta_rhodz=f_theta_rhodz(ind)
	454	theta_rhodzm1=f_theta_rhodzm1(ind)
	455	dtheta_rhodz=f_dtheta_rhodz(ind)
[129]	456
	457	IF (stage==1) THEN ! first stage : save model state in XXm1
[159]	458	DO l=ll_begin,ll_end
[174]	459	!$SIMD
	460	DO ij=ij_begin,ij_end
[148]	461	um1(ij+u_right,l)=u(ij+u_right,l)
	462	um1(ij+u_lup,l)=u(ij+u_lup,l)
	463	um1(ij+u_ldown,l)=u(ij+u_ldown,l)
	464	theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
	465	ENDDO
	466	ENDDO
[162]	467	END IF
[148]	468
[151]	469	DO l=ll_begin,ll_end
[174]	470	!$SIMD
	471	DO ij=ij_begin,ij_end
[148]	472	u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l)
	473	u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
	474	u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
	475	theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
	476	ENDDO
	477	ENDDO
[162]	478
[133]	479	IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
	480	hflux=f_hflux(ind); hfluxt=f_hfluxt(ind)
[138]	481	wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
	482	CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
[133]	483	END IF
[120]	484	END DO
[145]	485
	486	CALL trace_end("RK_scheme")
	487
[120]	488	END SUBROUTINE RK_scheme
	489
[129]	490	SUBROUTINE leapfrog_matsuno_scheme(stage)
[12]	491	IMPLICIT NONE
[129]	492	INTEGER :: ind, stage
	493	REAL :: tau
[145]	494
	495	CALL trace_start("leapfrog_matsuno_scheme")
	496
	497	tau = dt/nb_stage
[12]	498	DO ind=1,ndomain
[186]	499	IF (.NOT. assigned_domain(ind)) CYCLE
[138]	500	CALL swap_dimensions(ind)
	501	CALL swap_geometry(ind)
	502
[12]	503	ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
	504	psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
	505	psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
	506	dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
	507
	508
[129]	509	IF (stage==1) THEN
[12]	510	psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
[129]	511	ps(:)=psm1(:)+tau*dps(:)
	512	u(:,:)=um1(:,:)+tau*du(:,:)
	513	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
[12]	514
[129]	515	ELSE IF (stage==2) THEN
[12]	516
[129]	517	ps(:)=psm1(:)+tau*dps(:)
	518	u(:,:)=um1(:,:)+tau*du(:,:)
	519	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
[12]	520
	521	psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
	522	psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)
	523
	524	ELSE
	525
[129]	526	ps(:)=psm2(:)+2taudps(:)
	527	u(:,:)=um2(:,:)+2taudu(:,:)
	528	theta_rhodz(:,:)=theta_rhodzm2(:,:)+2taudtheta_rhodz(:,:)
[12]	529
	530	psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
	531	psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)
	532
	533	ENDIF
	534
	535	ENDDO
[145]	536	CALL trace_end("leapfrog_matsuno_scheme")
[12]	537
	538	END SUBROUTINE leapfrog_matsuno_scheme
	539
[50]	540	END SUBROUTINE timeloop
[133]	541
[159]	542	SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
[133]	543	USE icosa
[159]	544	USE omp_para
[133]	545	USE disvert_mod
[159]	546	IMPLICIT NONE
[133]	547	REAL(rstd), INTENT(IN) :: hflux(3iimjjm,llm), wflux(iim*jjm,llm+1)
	548	REAL(rstd), INTENT(INOUT) :: hfluxt(3iimjjm,llm), wfluxt(iim*jjm,llm+1)
	549	REAL(rstd), INTENT(IN) :: tau
[134]	550	LOGICAL, INTENT(INOUT) :: fluxt_zero
[148]	551	INTEGER :: l,i,j,ij
	552
[134]	553	IF(fluxt_zero) THEN
[151]	554
[134]	555	fluxt_zero=.FALSE.
[151]	556
	557	DO l=ll_begin,ll_end
[174]	558	!$SIMD
	559	DO ij=ij_begin_ext,ij_end_ext
[148]	560	hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l)
	561	hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l)
	562	hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l)
	563	ENDDO
	564	ENDDO
	565
[159]	566	IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
	567	DO l=ll_begin,ll_endp1
[174]	568	!$SIMD
	569	DO ij=ij_begin,ij_end
	570	wfluxt(ij,l) = tau*wflux(ij,l)
[159]	571	ENDDO
	572	ENDDO
	573	END IF
[162]	574
[134]	575	ELSE
[151]	576
	577	DO l=ll_begin,ll_end
[174]	578	!$SIMD
	579	DO ij=ij_begin_ext,ij_end_ext
[148]	580	hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l)
	581	hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l)
	582	hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l)
	583	ENDDO
	584	ENDDO
	585
[159]	586	IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
	587	DO l=ll_begin,ll_endp1
[174]	588	!$SIMD
	589	DO ij=ij_begin,ij_end
[159]	590	wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
	591	ENDDO
	592	ENDDO
	593	END IF
[148]	594
[159]	595	END IF
[151]	596
[133]	597	END SUBROUTINE accumulate_fluxes
[12]	598
[174]	599	! FUNCTION maxval_i(p)
	600	! USE icosa
	601	! IMPLICIT NONE
	602	! REAL(rstd), DIMENSION(iim*jjm) :: p
	603	! REAL(rstd) :: maxval_i
	604	! INTEGER :: j, ij
	605	!
	606	! maxval_i=p((jj_begin-1)*iim+ii_begin)
	607	!
	608	! DO j=jj_begin-1,jj_end+1
	609	! ij=(j-1)*iim
	610	! maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end)))
	611	! END DO
	612	! END FUNCTION maxval_i
[162]	613
[174]	614	! FUNCTION maxval_ik(p)
	615	! USE icosa
	616	! IMPLICIT NONE
	617	! REAL(rstd) :: p(iim*jjm, llm)
	618	! REAL(rstd) :: maxval_ik(llm)
	619	! INTEGER :: l,j, ij
	620	!
	621	! DO l=1,llm
	622	! maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l)
	623	! DO j=jj_begin-1,jj_end+1
	624	! ij=(j-1)*iim
	625	! maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l)))
	626	! END DO
	627	! END DO
	628	! END FUNCTION maxval_ik
[162]	629
[12]	630	END MODULE timeloop_gcm_mod

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

Download in other formats: