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

Last change on this file since 295 was 295, checked in by ymipsl, 10 years ago

Merging OpenMP parallisme mode : by subdomain and on vertical level.
This feature is actually experimental but may be retro-compatible with the last method based only on subdomain

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