Context Navigation

timeloop_gcm.f90 @ 3845

Last change on this file since 3845 was 3844, checked in by ymipsl, 10 years ago

Bug fix correction for some openMP configuration

File size: 19.0 KB

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

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

Download in other formats:

Original Format