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

Last change on this file since 173 was 171, checked in by ymipsl, 11 years ago
  • XIOS integration -

Compiling with "-with_xios" option. Adapt path to find XIOS library (arch.path)
Retro-compatible with the old output. If xios is not present, dynamico will use the standard writefield function.
Need to have the iodef.xml configuration file in the exec directory

YM
File size: 19.9 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  USE output_field_mod
43  USE write_field
44  IMPLICIT NONE
45
46    CHARACTER(len=255) :: def
47
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
75   
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
101   IF (xios_output) itau_out=1
102
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
123    CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
124
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
130!       f_massm1 => f_mass
131    ELSE ! eta = Lagrangian vertical coordinate
132       CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
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
138
139    def='runge_kutta'
140    CALL getin('scheme',def)
141
142    SELECT CASE (TRIM(def))
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
154     ! Model state 2 time steps ago (MLF)
155        CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
156        CALL allocate_field(f_um2,field_u,type_real,llm)
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
167    CASE default
168       PRINT*,'Bad selector for variable scheme : <', TRIM(def),             &
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
180    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
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)
187    CALL init_message(f_mass,req_i0,req_mass0)
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
193 
194  SUBROUTINE timeloop
195  USE icosa
196  USE dissip_gcm_mod
197  USE disvert_mod
198  USE caldyn_mod
199  USE caldyn_gcm_mod, ONLY : req_ps, req_mass
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
209  USE xios_mod
210  USE output_field_mod
211  IMPLICIT NONE 
212    REAL(rstd),POINTER :: q(:,:,:)
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(:,:)
217    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
218
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.
223
224    CALL caldyn_BC(f_phis, f_wflux) ! set constant values in first/last interfaces
225
226!$OMP BARRIER
227  DO ind=1,ndomain
228     CALL swap_dimensions(ind)
229     CALL swap_geometry(ind)
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
236  END DO
237  fluxt_zero=.TRUE.
238
239  DO it=0,itaumax
240   
241    CALL xios_update_calendar(it)
242    IF (MOD(it,itau_sync)==0) THEN
243      CALL send_message(f_ps,req_ps0)
244      CALL send_message(f_mass,req_mass0)
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)
249      CALL wait_message(req_mass0)
250      CALL wait_message(req_theta_rhodz0) 
251      CALL wait_message(req_u0)
252      CALL wait_message(req_q0) 
253    ENDIF
254   
255!    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
256    IF (mod(it,itau_out)==0 ) THEN
257      CALL update_time_counter(dt*it)
258      CALL output_field("q",f_q)
259      CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
260    ENDIF
261   
262    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
263
264    DO stage=1,nb_stage
265       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
266            f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, &
267            f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
268       SELECT CASE (scheme)
269       CASE(euler)
270          CALL euler_scheme(.TRUE.)
271       CASE (rk4)
272          CALL rk_scheme(stage)
273       CASE (mlf)
274          CALL  leapfrog_matsuno_scheme(stage)
275       CASE DEFAULT
276          STOP
277       END SELECT
278    END DO
279
280    IF (MOD(it+1,itau_dissip)==0) THEN
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
292
293    IF(MOD(it+1,itau_adv)==0) THEN
294
295       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
296       fluxt_zero=.TRUE.
297
298       ! FIXME : check that rhodz is consistent with ps
299       IF (check) THEN
300         DO ind=1,ndomain
301            CALL swap_dimensions(ind)
302            CALL swap_geometry(ind)
303            rhodz=f_rhodz(ind); ps=f_ps(ind);
304            CALL compute_rhodz(.FALSE., ps, rhodz)   
305         END DO
306       ENDIF
307
308    END IF
309
310
311
312!----------------------------------------------------
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)
317    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
318    ENDDO
319
320 
321  CONTAINS
322
323    SUBROUTINE Euler_scheme(with_dps)
324    IMPLICIT NONE
325    LOGICAL :: with_dps
326    INTEGER :: ind
327    INTEGER :: i,j,ij,l
328    CALL trace_start("Euler_scheme") 
329
330    DO ind=1,ndomain
331       CALL swap_dimensions(ind)
332       CALL swap_geometry(ind)
333
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
338                DO j=jj_begin,jj_end
339                   DO i=ii_begin,ii_end
340                      ij=(j-1)*iim+i
341                      ps(ij)=ps(ij)+dt*dps(ij)
342                   ENDDO
343                ENDDO
344             ENDIF
345          ELSE ! update mass
346             mass=f_mass(ind) ; dmass=f_dmass(ind) ;             
347             DO l=1,llm
348                DO j=jj_begin,jj_end
349                   DO i=ii_begin,ii_end
350                      ij=(j-1)*iim+i
351                      mass(ij,l)=mass(ij,l)+dt*dmass(ij,l)
352                   ENDDO
353                ENDDO
354             END DO
355          END IF
356
357          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
358          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
359          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
360       END IF ! update ps/mass
361       
362       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
363       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
364
365       DO l=ll_begin,ll_end
366         DO j=jj_begin,jj_end
367           DO i=ii_begin,ii_end
368             ij=(j-1)*iim+i
369             u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l)
370             u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
371             u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
372             theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l)
373           ENDDO
374         ENDDO
375       ENDDO
376    ENDDO
377
378    CALL trace_end("Euler_scheme") 
379
380    END SUBROUTINE Euler_scheme
381
382    SUBROUTINE RK_scheme(stage)
383      IMPLICIT NONE
384      INTEGER :: ind, stage
385      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
386      REAL(rstd) :: tau
387      INTEGER :: i,j,ij,l
388 
389      CALL trace_start("RK_scheme") 
390
391      tau = dt*coef(stage)
392
393      ! if mass coordinate, deal with ps first on one core
394      IF(caldyn_eta==eta_mass) THEN
395         IF (omp_first) THEN
396
397            DO ind=1,ndomain
398               CALL swap_dimensions(ind)
399               CALL swap_geometry(ind)
400               ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind) 
401               
402               IF (stage==1) THEN ! first stage : save model state in XXm1
403                  DO j=jj_begin,jj_end
404                     DO i=ii_begin,ii_end
405                        ij=(j-1)*iim+i
406                        psm1(ij)=ps(ij)
407                     ENDDO
408                  ENDDO
409               ENDIF
410           
411               ! updates are of the form x1 := x0 + tau*f(x1)
412               DO j=jj_begin,jj_end
413                  DO i=ii_begin,ii_end
414                     ij=(j-1)*iim+i
415                     ps(ij)=psm1(ij)+tau*dps(ij)
416                  ENDDO
417               ENDDO
418            ENDDO
419         ENDIF
420         CALL send_message(f_ps,req_ps)
421     
422      ELSE ! Lagrangian coordinate, deal with mass
423         DO ind=1,ndomain
424            CALL swap_dimensions(ind)
425            CALL swap_geometry(ind)
426            mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind)
427
428            IF (stage==1) THEN ! first stage : save model state in XXm1
429               DO l=ll_begin,ll_end
430                  DO j=jj_begin,jj_end
431                     DO i=ii_begin,ii_end
432                        ij=(j-1)*iim+i
433                        massm1(ij,l)=mass(ij,l)
434                     ENDDO
435                  ENDDO
436               ENDDO
437            END IF
438
439            ! updates are of the form x1 := x0 + tau*f(x1)
440            DO l=ll_begin,ll_end
441               DO j=jj_begin,jj_end
442                  DO i=ii_begin,ii_end
443                     ij=(j-1)*iim+i
444                     mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
445                  ENDDO
446               ENDDO
447            ENDDO
448         END DO
449         CALL send_message(f_mass,req_mass)
450
451      END IF
452
453      ! now deal with other prognostic variables
454      DO ind=1,ndomain
455         CALL swap_dimensions(ind)
456         CALL swap_geometry(ind)
457         u=f_u(ind)      ; du=f_du(ind)      ; um1=f_um1(ind) 
458         theta_rhodz=f_theta_rhodz(ind)
459         theta_rhodzm1=f_theta_rhodzm1(ind)
460         dtheta_rhodz=f_dtheta_rhodz(ind)
461         
462         IF (stage==1) THEN ! first stage : save model state in XXm1
463           DO l=ll_begin,ll_end
464             DO j=jj_begin,jj_end
465               DO i=ii_begin,ii_end
466                 ij=(j-1)*iim+i
467                 um1(ij+u_right,l)=u(ij+u_right,l)
468                 um1(ij+u_lup,l)=u(ij+u_lup,l)
469                 um1(ij+u_ldown,l)=u(ij+u_ldown,l)
470                 theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
471               ENDDO
472             ENDDO
473           ENDDO
474         END IF       
475
476         DO l=ll_begin,ll_end
477           DO j=jj_begin,jj_end
478             DO i=ii_begin,ii_end
479               ij=(j-1)*iim+i
480               u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l)
481               u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
482               u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
483               theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
484             ENDDO
485           ENDDO
486         ENDDO
487
488         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
489            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
490            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
491            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
492         END IF
493      END DO
494     
495      CALL trace_end("RK_scheme")
496     
497    END SUBROUTINE RK_scheme
498
499    SUBROUTINE leapfrog_matsuno_scheme(stage)
500    IMPLICIT NONE
501    INTEGER :: ind, stage
502    REAL :: tau
503
504      CALL trace_start("leapfrog_matsuno_scheme")
505   
506      tau = dt/nb_stage
507      DO ind=1,ndomain
508        CALL swap_dimensions(ind)
509        CALL swap_geometry(ind)
510
511        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
512        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
513        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
514        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
515
516       
517        IF (stage==1) THEN
518          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
519          ps(:)=psm1(:)+tau*dps(:)
520          u(:,:)=um1(:,:)+tau*du(:,:)
521          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
522
523        ELSE IF (stage==2) THEN
524
525          ps(:)=psm1(:)+tau*dps(:)
526          u(:,:)=um1(:,:)+tau*du(:,:)
527          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
528
529          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
530          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
531
532        ELSE
533
534          ps(:)=psm2(:)+2*tau*dps(:)
535          u(:,:)=um2(:,:)+2*tau*du(:,:)
536          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
537
538          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
539          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
540
541        ENDIF
542     
543      ENDDO
544      CALL trace_end("leapfrog_matsuno_scheme")
545     
546    END SUBROUTINE leapfrog_matsuno_scheme 
547         
548  END SUBROUTINE timeloop   
549
550 SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
551    USE icosa
552    USE omp_para
553    USE disvert_mod
554    IMPLICIT NONE
555    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
556    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
557    REAL(rstd), INTENT(IN) :: tau
558    LOGICAL, INTENT(INOUT) :: fluxt_zero
559    INTEGER :: l,i,j,ij
560
561    IF(fluxt_zero) THEN
562
563       fluxt_zero=.FALSE.
564
565       DO l=ll_begin,ll_end
566         DO j=jj_begin-1,jj_end+1
567           DO i=ii_begin-1,ii_end+1
568             ij=(j-1)*iim+i
569             hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l)
570             hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l)
571             hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l)
572           ENDDO
573         ENDDO
574       ENDDO
575
576       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
577          DO l=ll_begin,ll_endp1
578             DO j=jj_begin,jj_end
579                DO i=ii_begin,ii_end
580                   ij=(j-1)*iim+i
581                   wfluxt(ij,l) = tau*wflux(ij,l)
582                ENDDO
583             ENDDO
584          ENDDO
585       END IF
586
587    ELSE
588
589       DO l=ll_begin,ll_end
590         DO j=jj_begin-1,jj_end+1
591           DO i=ii_begin-1,ii_end+1
592             ij=(j-1)*iim+i
593             hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l)
594             hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l)
595             hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l)
596           ENDDO
597         ENDDO
598       ENDDO
599
600       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
601          DO l=ll_begin,ll_endp1
602             DO j=jj_begin,jj_end
603                DO i=ii_begin,ii_end
604                   ij=(j-1)*iim+i
605                   wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
606                ENDDO
607             ENDDO
608          ENDDO
609       END IF
610
611   END IF
612
613  END SUBROUTINE accumulate_fluxes
614 
615  FUNCTION maxval_i(p)
616    USE icosa
617    IMPLICIT NONE
618    REAL(rstd), DIMENSION(iim*jjm) :: p
619    REAL(rstd) :: maxval_i
620    INTEGER :: j, ij
621   
622    maxval_i=p((jj_begin-1)*iim+ii_begin)
623   
624    DO j=jj_begin-1,jj_end+1
625       ij=(j-1)*iim
626       maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end)))
627    END DO
628  END FUNCTION maxval_i
629
630  FUNCTION maxval_ik(p)
631    USE icosa
632    IMPLICIT NONE
633    REAL(rstd) :: p(iim*jjm, llm)
634    REAL(rstd) :: maxval_ik(llm)
635    INTEGER :: l,j, ij
636   
637    DO l=1,llm
638       maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l)
639       DO j=jj_begin-1,jj_end+1
640          ij=(j-1)*iim
641          maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l)))
642       END DO
643    END DO
644  END FUNCTION maxval_ik
645
646END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.