Context Navigation

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

Last change on this file since 276 was 266, checked in by ymipsl, 10 years ago

Synchronize trunk and Saturn branch.
Merge modification from Saturn branch to trunk

YM

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

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

Download in other formats: