Context Navigation

timeloop_gcm.f90 @ 314

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

activate splitting of XIOS file in physics so starting time is passed to the physic initialiszation.
call restart file periodically using the - itau_write_etat0 - start parameter.

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

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

Context Navigation

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

Download in other formats: