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

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

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