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

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