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

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

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