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

Last change on this file since 354 was 354, checked in by dubos, 9 years ago

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