Context Navigation

timeloop_gcm.f90 @ 3887

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

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

Download in other formats:

Original Format