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

Last change on this file since 178 was 174, checked in by ymipsl, 11 years ago

Transform 2 loops on i and j in one loop ij for efficiency, vectorization and future GPU programing

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