Context Navigation

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

Last change on this file since 195 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM

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

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

Download in other formats: