Context Navigation

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

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

Add the possibility to run without caldyn

=> scheme='none'

YM

File size: 19.4 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	ENDDO
306
307	CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q)
308
309	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
310
311	!$OMP MASTER
312	CALL SYSTEM_CLOCK(stop_clock)
313	CALL SYSTEM_CLOCK(count_rate=rate_clock)
314
315	IF (mpi_rank==0) THEN
316	PRINT ,"Time elapsed : ",(stop_clock-start_clock)1./rate_clock
317	ENDIF
318	!$OMP END MASTER
319
320	CONTAINS
321
322	SUBROUTINE Euler_scheme(with_dps)
323	IMPLICIT NONE
324	LOGICAL :: with_dps
325	INTEGER :: ind
326	INTEGER :: i,j,ij,l
327	CALL trace_start("Euler_scheme")
328
329	DO ind=1,ndomain
330	IF (.NOT. assigned_domain(ind)) CYCLE
331	CALL swap_dimensions(ind)
332	CALL swap_geometry(ind)
333
334	IF(with_dps) THEN ! update ps/mass
335	IF(caldyn_eta==eta_mass) THEN ! update ps
336	ps=f_ps(ind) ; dps=f_dps(ind) ;
337	IF (omp_first) THEN
338	!$SIMD
339	DO ij=ij_begin,ij_end
340	ps(ij)=ps(ij)+dt*dps(ij)
341	ENDDO
342	ENDIF
343	ELSE ! update mass
344	mass=f_mass(ind) ; dmass=f_dmass(ind) ;
345	DO l=1,llm
346	!$SIMD
347	DO ij=ij_begin,ij_end
348	mass(ij,l)=mass(ij,l)+dt*dmass(ij,l)
349	ENDDO
350	END DO
351	END IF
352
353	hflux=f_hflux(ind); hfluxt=f_hfluxt(ind)
354	wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
355	CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
356	END IF ! update ps/mass
357
358	u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
359	du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
360
361	DO l=ll_begin,ll_end
362	!$SIMD
363	DO ij=ij_begin,ij_end
364	u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l)
365	u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
366	u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
367	theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l)
368	ENDDO
369	ENDDO
370	ENDDO
371
372	CALL trace_end("Euler_scheme")
373
374	END SUBROUTINE Euler_scheme
375
376	SUBROUTINE RK_scheme(stage)
377	IMPLICIT NONE
378	INTEGER :: ind, stage
379	REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
380	REAL(rstd) :: tau
381	INTEGER :: i,j,ij,l
382
383	CALL trace_start("RK_scheme")
384
385	tau = dt*coef(stage)
386
387	! if mass coordinate, deal with ps first on one core
388	IF(caldyn_eta==eta_mass) THEN
389	IF (omp_first) THEN
390
391	DO ind=1,ndomain
392	IF (.NOT. assigned_domain(ind)) CYCLE
393	CALL swap_dimensions(ind)
394	CALL swap_geometry(ind)
395	ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind)
396
397	IF (stage==1) THEN ! first stage : save model state in XXm1
398	!$SIMD
399	DO ij=ij_begin,ij_end
400	psm1(ij)=ps(ij)
401	ENDDO
402	ENDIF
403
404	! updates are of the form x1 := x0 + tau*f(x1)
405	!$SIMD
406	DO ij=ij_begin,ij_end
407	ps(ij)=psm1(ij)+tau*dps(ij)
408	ENDDO
409	ENDDO
410	ENDIF
411	! CALL send_message(f_ps,req_ps)
412	!ym no overlap for now
413	! CALL wait_message(req_ps)
414
415	ELSE ! Lagrangian coordinate, deal with mass
416	DO ind=1,ndomain
417	IF (.NOT. assigned_domain(ind)) CYCLE
418	CALL swap_dimensions(ind)
419	CALL swap_geometry(ind)
420	mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind)
421
422	IF (stage==1) THEN ! first stage : save model state in XXm1
423	DO l=ll_begin,ll_end
424	!$SIMD
425	DO ij=ij_begin,ij_end
426	massm1(ij,l)=mass(ij,l)
427	ENDDO
428	ENDDO
429	END IF
430
431	! updates are of the form x1 := x0 + tau*f(x1)
432	DO l=ll_begin,ll_end
433	!$SIMD
434	DO ij=ij_begin,ij_end
435	mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
436	ENDDO
437	ENDDO
438	END DO
439	! CALL send_message(f_mass,req_mass)
440	!ym no overlap for now
441	! CALL wait_message(req_mass)
442
443	END IF
444
445	! now deal with other prognostic variables
446	DO ind=1,ndomain
447	IF (.NOT. assigned_domain(ind)) CYCLE
448	CALL swap_dimensions(ind)
449	CALL swap_geometry(ind)
450	u=f_u(ind) ; du=f_du(ind) ; um1=f_um1(ind)
451	theta_rhodz=f_theta_rhodz(ind)
452	theta_rhodzm1=f_theta_rhodzm1(ind)
453	dtheta_rhodz=f_dtheta_rhodz(ind)
454
455	IF (stage==1) THEN ! first stage : save model state in XXm1
456	DO l=ll_begin,ll_end
457	!$SIMD
458	DO ij=ij_begin,ij_end
459	um1(ij+u_right,l)=u(ij+u_right,l)
460	um1(ij+u_lup,l)=u(ij+u_lup,l)
461	um1(ij+u_ldown,l)=u(ij+u_ldown,l)
462	theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
463	ENDDO
464	ENDDO
465	END IF
466
467	DO l=ll_begin,ll_end
468	!$SIMD
469	DO ij=ij_begin,ij_end
470	u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l)
471	u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
472	u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
473	theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
474	ENDDO
475	ENDDO
476
477	IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
478	hflux=f_hflux(ind); hfluxt=f_hfluxt(ind)
479	wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
480	CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
481	END IF
482	END DO
483
484	CALL trace_end("RK_scheme")
485
486	END SUBROUTINE RK_scheme
487
488	SUBROUTINE leapfrog_matsuno_scheme(stage)
489	IMPLICIT NONE
490	INTEGER :: ind, stage
491	REAL :: tau
492
493	CALL trace_start("leapfrog_matsuno_scheme")
494
495	tau = dt/nb_stage
496	DO ind=1,ndomain
497	IF (.NOT. assigned_domain(ind)) CYCLE
498	CALL swap_dimensions(ind)
499	CALL swap_geometry(ind)
500
501	ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
502	psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
503	psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
504	dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
505
506
507	IF (stage==1) THEN
508	psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
509	ps(:)=psm1(:)+tau*dps(:)
510	u(:,:)=um1(:,:)+tau*du(:,:)
511	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
512
513	ELSE IF (stage==2) THEN
514
515	ps(:)=psm1(:)+tau*dps(:)
516	u(:,:)=um1(:,:)+tau*du(:,:)
517	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
518
519	psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
520	psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)
521
522	ELSE
523
524	ps(:)=psm2(:)+2taudps(:)
525	u(:,:)=um2(:,:)+2taudu(:,:)
526	theta_rhodz(:,:)=theta_rhodzm2(:,:)+2taudtheta_rhodz(:,:)
527
528	psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
529	psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)
530
531	ENDIF
532
533	ENDDO
534	CALL trace_end("leapfrog_matsuno_scheme")
535
536	END SUBROUTINE leapfrog_matsuno_scheme
537
538	END SUBROUTINE timeloop
539
540	SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
541	USE icosa
542	USE omp_para
543	USE disvert_mod
544	IMPLICIT NONE
545	REAL(rstd), INTENT(IN) :: hflux(3iimjjm,llm), wflux(iim*jjm,llm+1)
546	REAL(rstd), INTENT(INOUT) :: hfluxt(3iimjjm,llm), wfluxt(iim*jjm,llm+1)
547	REAL(rstd), INTENT(IN) :: tau
548	LOGICAL, INTENT(INOUT) :: fluxt_zero
549	INTEGER :: l,i,j,ij
550
551	IF(fluxt_zero) THEN
552
553	fluxt_zero=.FALSE.
554
555	DO l=ll_begin,ll_end
556	!$SIMD
557	DO ij=ij_begin_ext,ij_end_ext
558	hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l)
559	hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l)
560	hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l)
561	ENDDO
562	ENDDO
563
564	IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
565	DO l=ll_begin,ll_endp1
566	!$SIMD
567	DO ij=ij_begin,ij_end
568	wfluxt(ij,l) = tau*wflux(ij,l)
569	ENDDO
570	ENDDO
571	END IF
572
573	ELSE
574
575	DO l=ll_begin,ll_end
576	!$SIMD
577	DO ij=ij_begin_ext,ij_end_ext
578	hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l)
579	hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l)
580	hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l)
581	ENDDO
582	ENDDO
583
584	IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
585	DO l=ll_begin,ll_endp1
586	!$SIMD
587	DO ij=ij_begin,ij_end
588	wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
589	ENDDO
590	ENDDO
591	END IF
592
593	END IF
594
595	END SUBROUTINE accumulate_fluxes
596
597	! FUNCTION maxval_i(p)
598	! USE icosa
599	! IMPLICIT NONE
600	! REAL(rstd), DIMENSION(iim*jjm) :: p
601	! REAL(rstd) :: maxval_i
602	! INTEGER :: j, ij
603	!
604	! maxval_i=p((jj_begin-1)*iim+ii_begin)
605	!
606	! DO j=jj_begin-1,jj_end+1
607	! ij=(j-1)*iim
608	! maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end)))
609	! END DO
610	! END FUNCTION maxval_i
611
612	! FUNCTION maxval_ik(p)
613	! USE icosa
614	! IMPLICIT NONE
615	! REAL(rstd) :: p(iim*jjm, llm)
616	! REAL(rstd) :: maxval_ik(llm)
617	! INTEGER :: l,j, ij
618	!
619	! DO l=1,llm
620	! maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l)
621	! DO j=jj_begin-1,jj_end+1
622	! ij=(j-1)*iim
623	! maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l)))
624	! END DO
625	! END DO
626	! END FUNCTION maxval_ik
627
628	END MODULE timeloop_gcm_mod

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

Download in other formats: