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

Last change on this file since 327 was 327, checked in by ymipsl, 9 years ago

Merge recent developments from saturn branch onto trunk.

  • lmdz generic physics interface
  • performance improvment on mix mpi/openmp
  • asynchrone and overlaping communication
  • best domain distribution between process and threads
  • ....

This version is compatible with the actual saturn version and the both branches are considered merged on dynamico component.

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