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

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

Synced with aquaplanet branch HEAT@45 - tested with DCMIP4
File size: 19.3 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, 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
25CONTAINS
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 :",dt*It
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(:)+2*tau*dps(:)
552          u(:,:)=um2(:,:)+2*tau*du(:,:)
553          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_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(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
573    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,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 
624END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.