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

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

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