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

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

Time-accumulation of mass fluxes active
File size: 12.8 KB
Line 
1MODULE timeloop_gcm_mod
2
3  PRIVATE
4 
5  PUBLIC :: timeloop
6
7  INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
8
9CONTAINS
10 
11  SUBROUTINE timeloop
12  USE icosa
13  USE dissip_gcm_mod
14  USE caldyn_mod
15  USE etat0_mod
16  USE guided_mod
17  USE advect_tracer_mod
18  USE physics_mod
19  USE mpipara
20 
21  IMPLICIT NONE
22  TYPE(t_field),POINTER :: f_phis(:)
23!  TYPE(t_field),POINTER :: f_theta(:)
24  TYPE(t_field),POINTER :: f_q(:)
25  TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:)
26  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
27  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
28  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
29  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
30  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
31  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
32  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)
33
34  REAL(rstd),POINTER :: phis(:)
35  REAL(rstd),POINTER :: q(:,:,:)
36  REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
37  REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
38  REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
39  REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
40  REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
41  REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
42  REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
43
44!  REAL(rstd) :: dt, run_length
45  INTEGER :: ind
46  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
47  CHARACTER(len=255) :: scheme_name
48  LOGICAL :: fluxt_zero ! set to .TRUE. to start accumulating fluxes in time
49!  INTEGER :: itaumax
50!  REAL(rstd) ::write_period
51!  INTEGER    :: itau_out
52 
53!  dt=90.
54!  CALL getin('dt',dt)
55 
56!  itaumax=100
57!  CALL getin('itaumax',itaumax)
58
59!  run_length=dt*itaumax
60!  CALL getin('run_length',run_length)
61!  itaumax=run_length/dt
62!  PRINT *,'itaumax=',itaumax
63!  dt=dt/scale_factor
64
65! Trends
66  CALL allocate_field(f_dps,field_t,type_real)
67  CALL allocate_field(f_du,field_u,type_real,llm)
68  CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
69! Model state at current time step (RK/MLF/Euler)
70  CALL allocate_field(f_phis,field_t,type_real)
71  CALL allocate_field(f_ps,field_t,type_real)
72  CALL allocate_field(f_u,field_u,type_real,llm)
73  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
74! Model state at previous time step (RK/MLF)
75  CALL allocate_field(f_psm1,field_t,type_real)
76  CALL allocate_field(f_um1,field_u,type_real,llm)
77  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
78! Tracers
79  CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
80  CALL allocate_field(f_rhodz,field_t,type_real,llm)
81! Mass fluxes
82  CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
83  CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
84  CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
85  CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
86
87  scheme_name='runge_kutta'
88  CALL getin('scheme',scheme_name)
89
90  SELECT CASE (TRIM(scheme_name))
91  CASE('euler')
92     scheme=euler
93     nb_stage=1
94  CASE ('runge_kutta')
95     scheme=rk4
96     nb_stage=4
97  CASE ('leapfrog_matsuno')
98     scheme=mlf
99     matsuno_period=5
100     CALL getin('matsuno_period',matsuno_period)
101     nb_stage=matsuno_period+1
102     ! Model state 2 time steps ago (MLF)
103     CALL allocate_field(f_psm2,field_t,type_real)
104     CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
105     CALL allocate_field(f_um2,field_u,type_real,llm)
106  CASE default
107     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
108          ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
109     STOP
110  END SELECT
111 
112!  write_period=0
113!  CALL getin('write_period',write_period)
114!  write_period=write_period/scale_factor
115!  itau_out=FLOOR(.5+write_period/dt)
116!  PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
117 
118! Trends at previous time steps needed only by Adams-Bashforth
119!  CALL allocate_field(f_dpsm1,field_t,type_real)
120!  CALL allocate_field(f_dpsm2,field_t,type_real)
121!  CALL allocate_field(f_dum1,field_u,type_real,llm)
122!  CALL allocate_field(f_dum2,field_u,type_real,llm)
123!  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
124!  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
125
126!  CALL allocate_field(f_theta,field_t,type_real,llm)
127!  CALL allocate_field(f_dtheta,field_t,type_real,llm)
128
129  CALL init_dissip
130  CALL init_caldyn
131  CALL init_guided
132  CALL init_advect_tracer
133  CALL init_physics
134 
135  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
136  CALL writefield("phis",f_phis,once=.TRUE.)
137  CALL transfert_request(f_q,req_i1) 
138
139  DO ind=1,ndomain
140     CALL swap_dimensions(ind)
141     CALL swap_geometry(ind)
142     rhodz=f_rhodz(ind); ps=f_ps(ind)
143     CALL compute_rhodz(ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
144  END DO
145  fluxt_zero=.FALSE.
146
147  DO it=0,itaumax
148
149    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
150    IF (mod(it,itau_out)==0 ) THEN
151      CALL writefield("q",f_q)
152      CALL update_time_counter(dt*it)
153    ENDIF
154   
155!    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
156
157    DO stage=1,nb_stage
158       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
159            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
160            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
161       SELECT CASE (scheme)
162       CASE(euler)
163          CALL euler_scheme(.TRUE.)
164       CASE (rk4)
165          CALL rk_scheme(stage)
166       CASE (mlf)
167          CALL  leapfrog_matsuno_scheme(stage)
168         
169          !      CASE ('leapfrog')
170          !        CALL leapfrog_scheme
171          !
172          !      CASE ('adam_bashforth')
173          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
174          !        CALL adam_bashforth_scheme
175       CASE DEFAULT
176          STOP
177       END SELECT
178    END DO
179
180    CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
181    CALL euler_scheme(.FALSE.)
182
183    IF(MOD(it+1,itau_adv)==0) THEN
184       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
185       fluxt_zero=.TRUE.
186    END IF
187
188!    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
189    ENDDO
190 
191  CONTAINS
192
193    SUBROUTINE Euler_scheme(with_dps)
194    IMPLICIT NONE
195    LOGICAL :: with_dps
196    INTEGER :: ind
197     
198    DO ind=1,ndomain
199       IF(with_dps) THEN
200          ps=f_ps(ind) ; dps=f_dps(ind) ; 
201          ps(:)=ps(:)+dt*dps(:)
202          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
203          wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
204          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero)
205       END IF
206       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
207       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
208       u(:,:)=u(:,:)+dt*du(:,:)
209       theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
210    ENDDO
211
212    END SUBROUTINE Euler_scheme
213
214    SUBROUTINE RK_scheme(stage)
215      IMPLICIT NONE
216      INTEGER :: ind, stage
217      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
218      REAL(rstd) :: tau
219
220      tau = dt*coef(stage)
221     
222      DO ind=1,ndomain
223         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
224         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
225         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
226         
227         IF (stage==1) THEN ! first stage : save model state in XXm1
228            psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
229         END IF
230         ! updates are of the form x1 := x0 + tau*f(x1)
231         ps(:)=psm1(:)+tau*dps(:)
232         u(:,:)=um1(:,:)+tau*du(:,:)
233         theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
234         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
235            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
236            wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
237            CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero)
238         END IF
239      END DO
240    END SUBROUTINE RK_scheme
241
242    SUBROUTINE leapfrog_scheme
243    IMPLICIT NONE
244      INTEGER :: ind
245
246      DO ind=1,ndomain
247        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
248        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
249        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
250        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
251
252        IF (it==0) THEN
253          psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) 
254
255          ps(:)=ps(:)+dt*dps(:)
256          u(:,:)=u(:,:)+dt*du(:,:)
257          theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
258
259          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
260        ELSE
261       
262          ps(:)=psm2(:)+2*dt*dps(:)
263          u(:,:)=um2(:,:)+2*dt*du(:,:)
264          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
265
266          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
267          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
268        ENDIF
269         
270      ENDDO
271    END SUBROUTINE leapfrog_scheme 
272 
273    SUBROUTINE leapfrog_matsuno_scheme(stage)
274    IMPLICIT NONE
275    INTEGER :: ind, stage
276    REAL :: tau
277    tau = dt/nb_stage
278      DO ind=1,ndomain
279        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
280        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
281        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
282        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
283
284       
285!        IF (MOD(it,matsuno_period+1)==0) THEN
286        IF (stage==1) THEN
287          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
288          ps(:)=psm1(:)+tau*dps(:)
289          u(:,:)=um1(:,:)+tau*du(:,:)
290          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
291
292!        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
293        ELSE IF (stage==2) THEN
294
295          ps(:)=psm1(:)+tau*dps(:)
296          u(:,:)=um1(:,:)+tau*du(:,:)
297          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
298
299          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
300          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
301
302        ELSE
303
304          ps(:)=psm2(:)+2*tau*dps(:)
305          u(:,:)=um2(:,:)+2*tau*du(:,:)
306          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
307
308          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
309          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
310
311        ENDIF
312     
313      ENDDO
314     
315    END SUBROUTINE leapfrog_matsuno_scheme 
316         
317    SUBROUTINE adam_bashforth_scheme
318    IMPLICIT NONE
319      INTEGER :: ind
320
321      DO ind=1,ndomain
322        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
323        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
324        dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind)
325        dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind)
326
327        IF (it==0) THEN
328          dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
329          dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
330        ENDIF
331             
332        ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12
333        u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12
334        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12
335
336        dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
337        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
338
339      ENDDO     
340     
341    END SUBROUTINE adam_bashforth_scheme
342
343  END SUBROUTINE timeloop   
344
345  SUBROUTINE compute_rhodz(ps, rhodz)
346    USE icosa
347    USE disvert_mod
348    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
349    REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm)
350    INTEGER :: l,i,j,ij
351    DO l = 1, llm
352       DO j=jj_begin-1,jj_end+1
353          DO i=ii_begin-1,ii_end+1
354             ij=(j-1)*iim+i
355             rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
356          ENDDO
357       ENDDO
358    ENDDO
359  END SUBROUTINE compute_rhodz
360
361  SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau,fluxt_zero)
362    USE icosa
363    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
364    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
365    REAL(rstd), INTENT(IN) :: tau
366    LOGICAL, INTENT(INOUT) :: fluxt_zero
367    IF(fluxt_zero) THEN
368       fluxt_zero=.FALSE.
369       hfluxt = tau*hflux
370       wfluxt = tau*wflux
371    ELSE
372       hfluxt = hfluxt + tau*hflux
373       wfluxt = wfluxt + tau*wflux
374    END IF
375  END SUBROUTINE accumulate_fluxes
376 
377END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.