Context Navigation

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

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

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

Download in other formats: