Context Navigation

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

Last change on this file since 215 was 202, checked in by dubos, 10 years ago
OpenMP fixes (not tested)
File size: 20.6 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	IF (.NOT. assigned_domain(ind)) CYCLE
235	CALL swap_dimensions(ind)
236	CALL swap_geometry(ind)
237	rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
238	IF(caldyn_eta==eta_mass) THEN
239	CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
240	ELSE
241	rhodz(:,:)=mass(:,:)
242	END IF
243	END DO
244	fluxt_zero=.TRUE.
245
246	!$OMP MASTER
247	CALL SYSTEM_CLOCK(start_clock)
248	!$OMP END MASTER
249
250	CALL trace_on
251
252	DO it=0,itaumax
253
254	IF (xios_output) CALL xios_update_calendar(it)
255	IF (MOD(it,itau_sync)==0) THEN
256	CALL send_message(f_ps,req_ps0)
257	CALL wait_message(req_ps0)
258	CALL send_message(f_mass,req_mass0)
259	CALL wait_message(req_mass0)
260	CALL send_message(f_theta_rhodz,req_theta_rhodz0)
261	CALL wait_message(req_theta_rhodz0)
262	CALL send_message(f_u,req_u0)
263	CALL wait_message(req_u0)
264	CALL send_message(f_q,req_q0)
265	CALL wait_message(req_q0)
266
267	! CALL wait_message(req_ps0)
268	! CALL wait_message(req_mass0)
269	! CALL wait_message(req_theta_rhodz0)
270	! CALL wait_message(req_u0)
271	! CALL wait_message(req_q0)
272	ENDIF
273
274	!$OMP MASTER
275	IF (is_mpi_root) PRINT ,"It No :",It," t :",dtIt
276	!$OMP END MASTER
277	IF (mod(it,itau_out)==0 ) THEN
278	CALL update_time_counter(dt*it)
279	CALL output_field("q",f_q)
280	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
281	ENDIF
282
283	CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
284
285	DO stage=1,nb_stage
286	CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
287	f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, &
288	f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
289	SELECT CASE (scheme)
290	CASE(euler)
291	CALL euler_scheme(.TRUE.)
292	CASE (rk4)
293	CALL rk_scheme(stage)
294	CASE (mlf)
295	CALL leapfrog_matsuno_scheme(stage)
296	CASE DEFAULT
297	STOP
298	END SELECT
299	END DO
300
301	IF (MOD(it+1,itau_dissip)==0) THEN
302	! CALL send_message(f_ps,req_ps)
303	! CALL wait_message(req_ps)
304
305	IF(caldyn_eta==eta_mass) THEN
306	DO ind=1,ndomain
307	IF (.NOT. assigned_domain(ind)) CYCLE
308	CALL swap_dimensions(ind)
309	CALL swap_geometry(ind)
310	mass=f_mass(ind); ps=f_ps(ind);
311	CALL compute_rhodz(.TRUE., ps, mass)
312	END DO
313	ENDIF
314	! CALL send_message(f_mass,req_mass)
315	! CALL wait_message(req_mass)
316	CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz)
317	! CALL send_message(f_mass,req_mass)
318	! CALL wait_message(req_mass)
319	CALL euler_scheme(.FALSE.) ! update only u, theta
320	END IF
321
322	IF(MOD(it+1,itau_adv)==0) THEN
323
324	CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step
325	fluxt_zero=.TRUE.
326
327	! FIXME : check that rhodz is consistent with ps
328	IF (check) THEN
329	DO ind=1,ndomain
330	IF (.NOT. assigned_domain(ind)) CYCLE
331	CALL swap_dimensions(ind)
332	CALL swap_geometry(ind)
333	rhodz=f_rhodz(ind); ps=f_ps(ind);
334	CALL compute_rhodz(.FALSE., ps, rhodz)
335	END DO
336	ENDIF
337
338	END IF
339
340
341
342	!----------------------------------------------------
343	! jD_cur = jD_ref + day_ini - day_ref + it/day_step
344	! jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
345	! jD_cur = jD_cur + int(jH_cur)
346	! jH_cur = jH_cur - int(jH_cur)
347	CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
348
349	ENDDO
350
351	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
352
353	!$OMP MASTER
354	CALL SYSTEM_CLOCK(stop_clock)
355	CALL SYSTEM_CLOCK(count_rate=rate_clock)
356
357	IF (mpi_rank==0) THEN
358	PRINT ,"Time elapsed : ",(stop_clock-start_clock)1./rate_clock
359	ENDIF
360	!$OMP END MASTER
361
362	CONTAINS
363
364	SUBROUTINE Euler_scheme(with_dps)
365	IMPLICIT NONE
366	LOGICAL :: with_dps
367	INTEGER :: ind
368	INTEGER :: i,j,ij,l
369	CALL trace_start("Euler_scheme")
370
371	DO ind=1,ndomain
372	IF (.NOT. assigned_domain(ind)) CYCLE
373	CALL swap_dimensions(ind)
374	CALL swap_geometry(ind)
375
376	IF(with_dps) THEN ! update ps/mass
377	IF(caldyn_eta==eta_mass) THEN ! update ps
378	ps=f_ps(ind) ; dps=f_dps(ind) ;
379	IF (omp_first) THEN
380	!$SIMD
381	DO ij=ij_begin,ij_end
382	ps(ij)=ps(ij)+dt*dps(ij)
383	ENDDO
384	ENDIF
385	ELSE ! update mass
386	mass=f_mass(ind) ; dmass=f_dmass(ind) ;
387	DO l=1,llm
388	!$SIMD
389	DO ij=ij_begin,ij_end
390	mass(ij,l)=mass(ij,l)+dt*dmass(ij,l)
391	ENDDO
392	END DO
393	END IF
394
395	hflux=f_hflux(ind); hfluxt=f_hfluxt(ind)
396	wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
397	CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
398	END IF ! update ps/mass
399
400	u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
401	du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
402
403	DO l=ll_begin,ll_end
404	!$SIMD
405	DO ij=ij_begin,ij_end
406	u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l)
407	u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
408	u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
409	theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l)
410	ENDDO
411	ENDDO
412	ENDDO
413
414	CALL trace_end("Euler_scheme")
415
416	END SUBROUTINE Euler_scheme
417
418	SUBROUTINE RK_scheme(stage)
419	IMPLICIT NONE
420	INTEGER :: ind, stage
421	REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
422	REAL(rstd) :: tau
423	INTEGER :: i,j,ij,l
424
425	CALL trace_start("RK_scheme")
426
427	tau = dt*coef(stage)
428
429	! if mass coordinate, deal with ps first on one core
430	IF(caldyn_eta==eta_mass) THEN
431	IF (omp_first) THEN
432
433	DO ind=1,ndomain
434	IF (.NOT. assigned_domain(ind)) CYCLE
435	CALL swap_dimensions(ind)
436	CALL swap_geometry(ind)
437	ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind)
438
439	IF (stage==1) THEN ! first stage : save model state in XXm1
440	!$SIMD
441	DO ij=ij_begin,ij_end
442	psm1(ij)=ps(ij)
443	ENDDO
444	ENDIF
445
446	! updates are of the form x1 := x0 + tau*f(x1)
447	!$SIMD
448	DO ij=ij_begin,ij_end
449	ps(ij)=psm1(ij)+tau*dps(ij)
450	ENDDO
451	ENDDO
452	ENDIF
453	! CALL send_message(f_ps,req_ps)
454	!ym no overlap for now
455	! CALL wait_message(req_ps)
456
457	ELSE ! Lagrangian coordinate, deal with mass
458	DO ind=1,ndomain
459	IF (.NOT. assigned_domain(ind)) CYCLE
460	CALL swap_dimensions(ind)
461	CALL swap_geometry(ind)
462	mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind)
463
464	IF (stage==1) THEN ! first stage : save model state in XXm1
465	DO l=ll_begin,ll_end
466	!$SIMD
467	DO ij=ij_begin,ij_end
468	massm1(ij,l)=mass(ij,l)
469	ENDDO
470	ENDDO
471	END IF
472
473	! updates are of the form x1 := x0 + tau*f(x1)
474	DO l=ll_begin,ll_end
475	!$SIMD
476	DO ij=ij_begin,ij_end
477	mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
478	ENDDO
479	ENDDO
480	END DO
481	! CALL send_message(f_mass,req_mass)
482	!ym no overlap for now
483	! CALL wait_message(req_mass)
484
485	END IF
486
487	! now deal with other prognostic variables
488	DO ind=1,ndomain
489	IF (.NOT. assigned_domain(ind)) CYCLE
490	CALL swap_dimensions(ind)
491	CALL swap_geometry(ind)
492	u=f_u(ind) ; du=f_du(ind) ; um1=f_um1(ind)
493	theta_rhodz=f_theta_rhodz(ind)
494	theta_rhodzm1=f_theta_rhodzm1(ind)
495	dtheta_rhodz=f_dtheta_rhodz(ind)
496
497	IF (stage==1) THEN ! first stage : save model state in XXm1
498	DO l=ll_begin,ll_end
499	!$SIMD
500	DO ij=ij_begin,ij_end
501	um1(ij+u_right,l)=u(ij+u_right,l)
502	um1(ij+u_lup,l)=u(ij+u_lup,l)
503	um1(ij+u_ldown,l)=u(ij+u_ldown,l)
504	theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
505	ENDDO
506	ENDDO
507	END IF
508
509	DO l=ll_begin,ll_end
510	!$SIMD
511	DO ij=ij_begin,ij_end
512	u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l)
513	u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
514	u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
515	theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
516	ENDDO
517	ENDDO
518
519	IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
520	hflux=f_hflux(ind); hfluxt=f_hfluxt(ind)
521	wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
522	CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
523	END IF
524	END DO
525
526	CALL trace_end("RK_scheme")
527
528	END SUBROUTINE RK_scheme
529
530	SUBROUTINE leapfrog_matsuno_scheme(stage)
531	IMPLICIT NONE
532	INTEGER :: ind, stage
533	REAL :: tau
534
535	CALL trace_start("leapfrog_matsuno_scheme")
536
537	tau = dt/nb_stage
538	DO ind=1,ndomain
539	IF (.NOT. assigned_domain(ind)) CYCLE
540	CALL swap_dimensions(ind)
541	CALL swap_geometry(ind)
542
543	ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
544	psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
545	psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
546	dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
547
548
549	IF (stage==1) THEN
550	psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
551	ps(:)=psm1(:)+tau*dps(:)
552	u(:,:)=um1(:,:)+tau*du(:,:)
553	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
554
555	ELSE IF (stage==2) THEN
556
557	ps(:)=psm1(:)+tau*dps(:)
558	u(:,:)=um1(:,:)+tau*du(:,:)
559	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
560
561	psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
562	psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)
563
564	ELSE
565
566	ps(:)=psm2(:)+2taudps(:)
567	u(:,:)=um2(:,:)+2taudu(:,:)
568	theta_rhodz(:,:)=theta_rhodzm2(:,:)+2taudtheta_rhodz(:,:)
569
570	psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
571	psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)
572
573	ENDIF
574
575	ENDDO
576	CALL trace_end("leapfrog_matsuno_scheme")
577
578	END SUBROUTINE leapfrog_matsuno_scheme
579
580	END SUBROUTINE timeloop
581
582	SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
583	USE icosa
584	USE omp_para
585	USE disvert_mod
586	IMPLICIT NONE
587	REAL(rstd), INTENT(IN) :: hflux(3iimjjm,llm), wflux(iim*jjm,llm+1)
588	REAL(rstd), INTENT(INOUT) :: hfluxt(3iimjjm,llm), wfluxt(iim*jjm,llm+1)
589	REAL(rstd), INTENT(IN) :: tau
590	LOGICAL, INTENT(INOUT) :: fluxt_zero
591	INTEGER :: l,i,j,ij
592
593	IF(fluxt_zero) THEN
594
595	fluxt_zero=.FALSE.
596
597	DO l=ll_begin,ll_end
598	!$SIMD
599	DO ij=ij_begin_ext,ij_end_ext
600	hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l)
601	hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l)
602	hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l)
603	ENDDO
604	ENDDO
605
606	IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
607	DO l=ll_begin,ll_endp1
608	!$SIMD
609	DO ij=ij_begin,ij_end
610	wfluxt(ij,l) = tau*wflux(ij,l)
611	ENDDO
612	ENDDO
613	END IF
614
615	ELSE
616
617	DO l=ll_begin,ll_end
618	!$SIMD
619	DO ij=ij_begin_ext,ij_end_ext
620	hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l)
621	hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l)
622	hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l)
623	ENDDO
624	ENDDO
625
626	IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
627	DO l=ll_begin,ll_endp1
628	!$SIMD
629	DO ij=ij_begin,ij_end
630	wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
631	ENDDO
632	ENDDO
633	END IF
634
635	END IF
636
637	END SUBROUTINE accumulate_fluxes
638
639	! FUNCTION maxval_i(p)
640	! USE icosa
641	! IMPLICIT NONE
642	! REAL(rstd), DIMENSION(iim*jjm) :: p
643	! REAL(rstd) :: maxval_i
644	! INTEGER :: j, ij
645	!
646	! maxval_i=p((jj_begin-1)*iim+ii_begin)
647	!
648	! DO j=jj_begin-1,jj_end+1
649	! ij=(j-1)*iim
650	! maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end)))
651	! END DO
652	! END FUNCTION maxval_i
653
654	! FUNCTION maxval_ik(p)
655	! USE icosa
656	! IMPLICIT NONE
657	! REAL(rstd) :: p(iim*jjm, llm)
658	! REAL(rstd) :: maxval_ik(llm)
659	! INTEGER :: l,j, ij
660	!
661	! DO l=1,llm
662	! maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l)
663	! DO j=jj_begin-1,jj_end+1
664	! ij=(j-1)*iim
665	! maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l)))
666	! END DO
667	! END DO
668	! END FUNCTION maxval_ik
669
670	END MODULE timeloop_gcm_mod

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

Download in other formats: