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

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

Almost nothing
File size: 19.5 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!$OMP BARRIER
203  DO ind=1,ndomain
204     CALL swap_dimensions(ind)
205     CALL swap_geometry(ind)
206     rhodz=f_rhodz(ind); ps=f_ps(ind)
207     CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
208  END DO
209  fluxt_zero=.TRUE.
210
211  DO it=0,itaumax
212    IF (MOD(it,itau_sync)==0) THEN
213      CALL send_message(f_ps,req_ps0)
214      CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
215      CALL send_message(f_u,req_u0)
216      CALL send_message(f_q,req_q0) 
217      CALL wait_message(req_ps0)
218      CALL wait_message(req_theta_rhodz0) 
219      CALL wait_message(req_u0)
220      CALL wait_message(req_q0) 
221    ENDIF
222   
223!    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
224    IF (mod(it,itau_out)==0 ) THEN
225      CALL writefield("q",f_q)
226      CALL update_time_counter(dt*it)
227      CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
228    ENDIF
229   
230    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
231
232    DO stage=1,nb_stage
233       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
234            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
235            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
236       SELECT CASE (scheme)
237       CASE(euler)
238          CALL euler_scheme(.TRUE.)
239       CASE (rk4)
240          CALL rk_scheme(stage)
241       CASE (mlf)
242          CALL  leapfrog_matsuno_scheme(stage)
243         
244          !      CASE ('leapfrog')
245          !        CALL leapfrog_scheme
246          !
247          !      CASE ('adam_bashforth')
248          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
249          !        CALL adam_bashforth_scheme
250       CASE DEFAULT
251          STOP
252       END SELECT
253    END DO
254
255    IF (MOD(it+1,itau_dissip)==0) THEN
256      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
257      CALL euler_scheme(.FALSE.)
258    ENDIF
259
260    IF(MOD(it+1,itau_adv)==0) THEN
261
262       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
263       fluxt_zero=.TRUE.
264
265       ! FIXME : check that rhodz is consistent with ps
266       IF (check) THEN
267         DO ind=1,ndomain
268            CALL swap_dimensions(ind)
269            CALL swap_geometry(ind)
270            rhodz=f_rhodz(ind); ps=f_ps(ind);
271            CALL compute_rhodz(.FALSE., ps, rhodz)   
272         END DO
273       ENDIF
274
275    END IF
276
277
278
279!----------------------------------------------------
280!    jD_cur = jD_ref + day_ini - day_ref + it/day_step
281!    jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
282!    jD_cur = jD_cur + int(jH_cur)
283!    jH_cur = jH_cur - int(jH_cur)
284!    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
285
286!    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
287    ENDDO
288
289 
290  CONTAINS
291
292    SUBROUTINE Euler_scheme(with_dps)
293    IMPLICIT NONE
294    LOGICAL :: with_dps
295    INTEGER :: ind
296    INTEGER :: i,j,ij,l
297    CALL trace_start("Euler_scheme") 
298
299    DO ind=1,ndomain
300       CALL swap_dimensions(ind)
301       CALL swap_geometry(ind)
302       IF(with_dps) THEN
303         ps=f_ps(ind) ; dps=f_dps(ind) ; 
304
305         IF (omp_first) THEN
306           DO j=jj_begin,jj_end
307             DO i=ii_begin,ii_end
308               ij=(j-1)*iim+i
309               ps(ij)=ps(ij)+dt*dps(ij)
310             ENDDO
311           ENDDO
312         ENDIF
313         
314         hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
315         wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
316         CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
317       END IF
318       
319       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
320       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
321
322       DO l=ll_begin,ll_end
323         DO j=jj_begin,jj_end
324           DO i=ii_begin,ii_end
325             ij=(j-1)*iim+i
326             u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l)
327             u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
328             u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
329             theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l)
330           ENDDO
331         ENDDO
332       ENDDO
333    ENDDO
334
335    CALL trace_end("Euler_scheme") 
336
337    END SUBROUTINE Euler_scheme
338
339    SUBROUTINE RK_scheme(stage)
340      USE caldyn_gcm_mod
341      IMPLICIT NONE
342      INTEGER :: ind, stage
343      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
344      REAL(rstd) :: tau
345      INTEGER :: i,j,ij,l
346 
347      CALL trace_start("RK_scheme") 
348
349      tau = dt*coef(stage)
350
351      DO ind=1,ndomain
352         CALL swap_dimensions(ind)
353         CALL swap_geometry(ind)
354         ps=f_ps(ind)   
355         psm1=f_psm1(ind) 
356         dps=f_dps(ind) 
357         
358         IF (stage==1) THEN ! first stage : save model state in XXm1
359           IF (omp_first) THEN
360             DO j=jj_begin,jj_end
361               DO i=ii_begin,ii_end
362                 ij=(j-1)*iim+i
363                 psm1(ij)=ps(ij)
364               ENDDO
365             ENDDO
366           ENDIF
367         END IF
368         
369         ! updates are of the form x1 := x0 + tau*f(x1)
370           IF (omp_first) THEN
371             DO j=jj_begin,jj_end
372               DO i=ii_begin,ii_end
373                 ij=(j-1)*iim+i
374                 ps(ij)=psm1(ij)+tau*dps(ij)
375               ENDDO
376             ENDDO
377           ENDIF
378       
379       ENDDO
380
381       CALL send_message(f_ps,req_ps)     
382
383
384     
385      DO ind=1,ndomain
386         CALL swap_dimensions(ind)
387         CALL swap_geometry(ind)
388         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
389         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
390         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
391         
392         IF (stage==1) THEN ! first stage : save model state in XXm1
393             
394         DO l=ll_begin,ll_end
395             DO j=jj_begin,jj_end
396               DO i=ii_begin,ii_end
397                 ij=(j-1)*iim+i
398                 um1(ij+u_right,l)=u(ij+u_right,l)
399                 um1(ij+u_lup,l)=u(ij+u_lup,l)
400                 um1(ij+u_ldown,l)=u(ij+u_ldown,l)
401                 theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
402               ENDDO
403             ENDDO
404           ENDDO
405
406         END IF
407         ! updates are of the form x1 := x0 + tau*f(x1)
408         
409         DO l=ll_begin,ll_end
410           DO j=jj_begin,jj_end
411             DO i=ii_begin,ii_end
412               ij=(j-1)*iim+i
413               u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l)
414               u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
415               u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
416               theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
417             ENDDO
418           ENDDO
419         ENDDO
420
421         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
422            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
423            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
424            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
425         END IF
426      END DO
427     
428      CALL trace_end("RK_scheme")
429     
430    END SUBROUTINE RK_scheme
431
432    SUBROUTINE leapfrog_scheme
433    IMPLICIT NONE
434      INTEGER :: ind
435
436      CALL trace_start("leapfrog_scheme")
437       
438      DO ind=1,ndomain
439        CALL swap_dimensions(ind)
440        CALL swap_geometry(ind)
441        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
442        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
443        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
444        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
445
446        IF (it==0) THEN
447          psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) 
448
449          ps(:)=ps(:)+dt*dps(:)
450          u(:,:)=u(:,:)+dt*du(:,:)
451          theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
452
453          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
454        ELSE
455       
456          ps(:)=psm2(:)+2*dt*dps(:)
457          u(:,:)=um2(:,:)+2*dt*du(:,:)
458          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
459
460          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
461          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
462        ENDIF
463         
464      ENDDO
465
466      CALL trace_end("leapfrog_scheme")
467
468    END SUBROUTINE leapfrog_scheme 
469 
470    SUBROUTINE leapfrog_matsuno_scheme(stage)
471    IMPLICIT NONE
472    INTEGER :: ind, stage
473    REAL :: tau
474
475      CALL trace_start("leapfrog_matsuno_scheme")
476   
477      tau = dt/nb_stage
478      DO ind=1,ndomain
479        CALL swap_dimensions(ind)
480        CALL swap_geometry(ind)
481
482        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
483        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
484        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
485        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
486
487       
488        IF (stage==1) THEN
489          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
490          ps(:)=psm1(:)+tau*dps(:)
491          u(:,:)=um1(:,:)+tau*du(:,:)
492          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
493
494        ELSE IF (stage==2) THEN
495
496          ps(:)=psm1(:)+tau*dps(:)
497          u(:,:)=um1(:,:)+tau*du(:,:)
498          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
499
500          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
501          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
502
503        ELSE
504
505          ps(:)=psm2(:)+2*tau*dps(:)
506          u(:,:)=um2(:,:)+2*tau*du(:,:)
507          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
508
509          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
510          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
511
512        ENDIF
513     
514      ENDDO
515      CALL trace_end("leapfrog_matsuno_scheme")
516     
517    END SUBROUTINE leapfrog_matsuno_scheme 
518         
519    SUBROUTINE adam_bashforth_scheme
520    IMPLICIT NONE
521      INTEGER :: ind
522
523      CALL trace_start("adam_bashforth_scheme")
524
525      DO ind=1,ndomain
526        CALL swap_dimensions(ind)
527        CALL swap_geometry(ind)
528        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
529        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
530        dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind)
531        dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind)
532
533        IF (it==0) THEN
534          dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
535          dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
536        ENDIF
537             
538        ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12
539        u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12
540        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12
541
542        dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
543        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
544
545      ENDDO     
546
547      CALL trace_end("adam_bashforth_scheme")
548     
549    END SUBROUTINE adam_bashforth_scheme
550
551  END SUBROUTINE timeloop   
552
553  SUBROUTINE compute_rhodz(comp, ps, rhodz)
554    USE icosa
555    USE disvert_mod
556    USE omp_para
557    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
558    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
559    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
560    REAL(rstd) :: m, err
561    INTEGER :: l,i,j,ij,dd
562    err=0.
563
564    dd=0
565
566    DO l = ll_begin, ll_end
567       DO j=jj_begin-dd,jj_end+dd
568          DO i=ii_begin-dd,ii_end+dd
569             ij=(j-1)*iim+i
570             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
571             IF(comp) THEN
572                rhodz(ij,l) = m
573             ELSE
574                err = MAX(err,abs(m-rhodz(ij,l)))
575             END IF
576          ENDDO
577       ENDDO
578    ENDDO
579
580    IF(.NOT. comp) THEN
581       IF(err>1e-10) THEN
582          PRINT *, 'Discrepancy between ps and rhodz detected', err
583          STOP
584       ELSE
585!          PRINT *, 'No discrepancy between ps and rhodz detected'
586       END IF
587    END IF
588
589  END SUBROUTINE compute_rhodz
590
591  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
592    USE icosa
593    USE omp_para
594    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
595    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
596    REAL(rstd), INTENT(IN) :: tau
597    LOGICAL, INTENT(INOUT) :: fluxt_zero
598    INTEGER :: l,i,j,ij
599
600    IF(fluxt_zero) THEN
601
602       fluxt_zero=.FALSE.
603
604       DO l=ll_begin,ll_end
605         DO j=jj_begin-1,jj_end+1
606           DO i=ii_begin-1,ii_end+1
607             ij=(j-1)*iim+i
608             hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l)
609             hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l)
610             hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l)
611           ENDDO
612         ENDDO
613       ENDDO
614
615       DO l=ll_begin,ll_endp1
616         DO j=jj_begin,jj_end
617           DO i=ii_begin,ii_end
618             ij=(j-1)*iim+i
619             wfluxt(ij,l) = tau*wflux(ij,l)
620           ENDDO
621         ENDDO
622       ENDDO
623    ELSE
624
625       DO l=ll_begin,ll_end
626         DO j=jj_begin-1,jj_end+1
627           DO i=ii_begin-1,ii_end+1
628             ij=(j-1)*iim+i
629             hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l)
630             hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l)
631             hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l)
632           ENDDO
633         ENDDO
634       ENDDO
635
636       DO l=ll_begin,ll_endp1
637         DO j=jj_begin,jj_end
638           DO i=ii_begin,ii_end
639             ij=(j-1)*iim+i
640             wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
641           ENDDO
642         ENDDO
643       ENDDO
644
645    END IF
646
647  END SUBROUTINE accumulate_fluxes
648 
649END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.