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
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  :: itau_sync=10
10
11  TYPE(t_message) :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0
12
13  TYPE(t_field),POINTER :: f_q(:)
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(:)
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
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 ioipsl
42  IMPLICIT NONE
43
44    CHARACTER(len=255) :: def
45
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
98
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
111
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
132    CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
133
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
139!       f_massm1 => f_mass
140    ELSE ! eta = Lagrangian vertical coordinate
141       CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
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
147
148    def='runge_kutta'
149    CALL getin('scheme',def)
150
151    SELECT CASE (TRIM(def))
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
163     ! Model state 2 time steps ago (MLF)
164        CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
165        CALL allocate_field(f_um2,field_u,type_real,llm)
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
176    CASE default
177       PRINT*,'Bad selector for variable scheme : <', TRIM(def),             &
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
189    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
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)
196    CALL init_message(f_mass,req_i0,req_mass0)
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
202 
203  SUBROUTINE timeloop
204  USE icosa
205  USE dissip_gcm_mod
206  USE disvert_mod
207  USE caldyn_mod
208  USE caldyn_gcm_mod, ONLY : req_ps, req_mass
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(:,:,:)
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(:,:)
224    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
225
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.
230
231    CALL caldyn_BC(f_phis, f_wflux) ! set constant values in first/last interfaces
232
233!$OMP BARRIER
234  DO ind=1,ndomain
235     CALL swap_dimensions(ind)
236     CALL swap_geometry(ind)
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
243  END DO
244  fluxt_zero=.TRUE.
245
246  DO it=0,itaumax
247    IF (MOD(it,itau_sync)==0) THEN
248      CALL send_message(f_ps,req_ps0)
249      CALL send_message(f_mass,req_mass0)
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)
254      CALL wait_message(req_mass0)
255      CALL wait_message(req_theta_rhodz0) 
256      CALL wait_message(req_u0)
257      CALL wait_message(req_q0) 
258    ENDIF
259   
260!    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
261    IF (mod(it,itau_out)==0 ) THEN
262      CALL writefield("q",f_q)
263      CALL update_time_counter(dt*it)
264      CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
265    ENDIF
266   
267    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
268
269    DO stage=1,nb_stage
270       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
271            f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, &
272            f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
273       SELECT CASE (scheme)
274       CASE(euler)
275          CALL euler_scheme(.TRUE.)
276       CASE (rk4)
277          CALL rk_scheme(stage)
278       CASE (mlf)
279          CALL  leapfrog_matsuno_scheme(stage)
280       CASE DEFAULT
281          STOP
282       END SELECT
283    END DO
284
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
289
290    IF(MOD(it+1,itau_adv)==0) THEN
291
292       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
293       fluxt_zero=.TRUE.
294
295       ! FIXME : check that rhodz is consistent with ps
296       IF (check) THEN
297         DO ind=1,ndomain
298            CALL swap_dimensions(ind)
299            CALL swap_geometry(ind)
300            rhodz=f_rhodz(ind); ps=f_ps(ind);
301            CALL compute_rhodz(.FALSE., ps, rhodz)   
302         END DO
303       ENDIF
304
305    END IF
306
307
308
309!----------------------------------------------------
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
316!    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
317    ENDDO
318
319 
320  CONTAINS
321
322    SUBROUTINE Euler_scheme(with_dps)
323    IMPLICIT NONE
324    LOGICAL :: with_dps
325    INTEGER :: ind
326    INTEGER :: i,j,ij,l
327    CALL trace_start("Euler_scheme") 
328
329    DO ind=1,ndomain
330       CALL swap_dimensions(ind)
331       CALL swap_geometry(ind)
332
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
360       
361       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
362       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
363
364       DO l=ll_begin,ll_end
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
375    ENDDO
376
377    CALL trace_end("Euler_scheme") 
378
379    END SUBROUTINE Euler_scheme
380
381    SUBROUTINE RK_scheme(stage)
382      IMPLICIT NONE
383      INTEGER :: ind, stage
384      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
385      REAL(rstd) :: tau
386      INTEGER :: i,j,ij,l
387 
388      CALL trace_start("RK_scheme") 
389
390      tau = dt*coef(stage)
391
392      ! if mass coordinate, deal with ps first on one core
393      IF(caldyn_eta==eta_mass) THEN
394         IF (omp_first) THEN
395
396            DO ind=1,ndomain
397               CALL swap_dimensions(ind)
398               CALL swap_geometry(ind)
399               ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind) 
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
416               ENDDO
417            ENDDO
418         ENDIF
419         CALL send_message(f_ps,req_ps)
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
450      END IF
451
452      ! now deal with other prognostic variables
453      DO ind=1,ndomain
454         CALL swap_dimensions(ind)
455         CALL swap_geometry(ind)
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)
460         
461         IF (stage==1) THEN ! first stage : save model state in XXm1
462           DO l=ll_begin,ll_end
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
473         END IF       
474
475         DO l=ll_begin,ll_end
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
486
487         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
488            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
489            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
490            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
491         END IF
492      END DO
493     
494      CALL trace_end("RK_scheme")
495     
496    END SUBROUTINE RK_scheme
497
498    SUBROUTINE leapfrog_matsuno_scheme(stage)
499    IMPLICIT NONE
500    INTEGER :: ind, stage
501    REAL :: tau
502
503      CALL trace_start("leapfrog_matsuno_scheme")
504   
505      tau = dt/nb_stage
506      DO ind=1,ndomain
507        CALL swap_dimensions(ind)
508        CALL swap_geometry(ind)
509
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       
516        IF (stage==1) THEN
517          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
518          ps(:)=psm1(:)+tau*dps(:)
519          u(:,:)=um1(:,:)+tau*du(:,:)
520          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
521
522        ELSE IF (stage==2) THEN
523
524          ps(:)=psm1(:)+tau*dps(:)
525          u(:,:)=um1(:,:)+tau*du(:,:)
526          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
527
528          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
529          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
530
531        ELSE
532
533          ps(:)=psm2(:)+2*tau*dps(:)
534          u(:,:)=um2(:,:)+2*tau*du(:,:)
535          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
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
543      CALL trace_end("leapfrog_matsuno_scheme")
544     
545    END SUBROUTINE leapfrog_matsuno_scheme 
546         
547  END SUBROUTINE timeloop   
548
549 SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
550    USE icosa
551    USE omp_para
552    USE disvert_mod
553    IMPLICIT NONE
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
557    LOGICAL, INTENT(INOUT) :: fluxt_zero
558    INTEGER :: l,i,j,ij
559
560    IF(fluxt_zero) THEN
561
562       fluxt_zero=.FALSE.
563
564       DO l=ll_begin,ll_end
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
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
585
586    ELSE
587
588       DO l=ll_begin,ll_end
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
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
609
610   END IF
611
612  END SUBROUTINE accumulate_fluxes
613 
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
645END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.