Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/timeloop_gcm.f90 @ 268

Last change on this file since 268 was 262, checked in by ymipsl, 10 years ago

Save starting iteration itau0 in dynamico restart file.
Timeloop now begin from itau0 which can be read from start file or set to 0 by default.

YM

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

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

Download in other formats: