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

Last change on this file since 162 was 162, checked in by dubos, 11 years ago

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