Context Navigation

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

Rev	Line
[12]	1	MODULE timeloop_gcm_mod
[133]	2
	3	PRIVATE
[12]	4
[133]	5	PUBLIC :: timeloop
	6
	7	INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
[148]	8	INTEGER :: itau_sync=10
[133]	9
[12]	10	CONTAINS
	11
	12	SUBROUTINE timeloop
[19]	13	USE icosa
[15]	14	USE dissip_gcm_mod
[17]	15	USE caldyn_mod
[12]	16	USE etat0_mod
[17]	17	USE guided_mod
	18	USE advect_tracer_mod
[81]	19	USE physics_mod
[131]	20	USE mpipara
[149]	21	USE IOIPSL
	22	USE maxicosa
	23	USE check_conserve_mod
[145]	24	USE trace
[148]	25	USE transfert_mod
[12]	26	IMPLICIT NONE
	27	TYPE(t_field),POINTER :: f_phis(:)
[129]	28	! TYPE(t_field),POINTER :: f_theta(:)
[17]	29	TYPE(t_field),POINTER :: f_q(:)
[132]	30	TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:)
[12]	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(:)
[133]	37	TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)
[12]	38
	39	REAL(rstd),POINTER :: phis(:)
[17]	40	REAL(rstd),POINTER :: q(:,:,:)
[12]	41	REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
	42	REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
[133]	43	REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
[12]	44	REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
	45	REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
	46	REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
[133]	47	REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
[98]	48	! REAL(rstd) :: dt, run_length
[12]	49	INTEGER :: ind
[133]	50	INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
	51	CHARACTER(len=255) :: scheme_name
[138]	52	LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
[149]	53	CHARACTER(len=7) :: first
	54	REAL(rstd),SAVE :: jD_cur, jH_cur
	55	REAL(rstd),SAVE :: start_time
[148]	56	LOGICAL, PARAMETER :: check=.FALSE.
[98]	57	! INTEGER :: itaumax
	58	! REAL(rstd) ::write_period
	59	! INTEGER :: itau_out
[63]	60
[98]	61	! dt=90.
	62	! CALL getin('dt',dt)
[32]	63
[98]	64	! itaumax=100
	65	! CALL getin('itaumax',itaumax)
[17]	66
[98]	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
[45]	72
[129]	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)
[132]	88	CALL allocate_field(f_rhodz,field_t,type_real,llm)
[134]	89	! Mass fluxes
[133]	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)
[129]	94
[149]	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
[133]	147	scheme_name='runge_kutta'
	148	CALL getin('scheme',scheme_name)
[129]	149
[133]	150	SELECT CASE (TRIM(scheme_name))
[129]	151	CASE('euler')
[133]	152	scheme=euler
[129]	153	nb_stage=1
	154	CASE ('runge_kutta')
[133]	155	scheme=rk4
[129]	156	nb_stage=4
	157	CASE ('leapfrog_matsuno')
[133]	158	scheme=mlf
[129]	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
[133]	167	PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name), &
[129]	168	' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
	169	STOP
	170	END SELECT
[12]	171
[98]	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
[63]	177
[129]	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)
[12]	187
[98]	188	CALL init_dissip
	189	CALL init_caldyn
	190	CALL init_guided
	191	CALL init_advect_tracer
	192	CALL init_physics
[149]	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)
[97]	196	CALL writefield("phis",f_phis,once=.TRUE.)
[81]	197	CALL transfert_request(f_q,req_i1)
[50]	198
[133]	199	DO ind=1,ndomain
	200	CALL swap_dimensions(ind)
	201	CALL swap_geometry(ind)
	202	rhodz=f_rhodz(ind); ps=f_ps(ind)
[138]	203	CALL compute_rhodz(.TRUE., ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
[133]	204	END DO
[138]	205	fluxt_zero=.TRUE.
[132]	206
[138]	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
[149]	216
[148]	217	CALL transfert_request(f_phis,req_i0)
	218	CALL transfert_request(f_phis,req_i1)
[149]	219	CALL transfert_request(f_phis,req_i1)
[148]	220
[12]	221	DO it=0,itaumax
[148]	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
[63]	229	IF (mod(it,itau_out)==0 ) THEN
[149]	230	! IF (is_mpi_root) PRINT ,"It No :",It," t :",dtIt
[81]	231	CALL update_time_counter(dt*it)
[149]	232	CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
[63]	233	ENDIF
[129]	234
[149]	235	CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
	236
[129]	237	DO stage=1,nb_stage
	238	CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
[132]	239	f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
[134]	240	f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
[133]	241	SELECT CASE (scheme)
	242	CASE(euler)
	243	CALL euler_scheme(.TRUE.)
	244	CASE (rk4)
[129]	245	CALL rk_scheme(stage)
[133]	246	CASE (mlf)
[129]	247	CALL leapfrog_matsuno_scheme(stage)
	248	! CASE ('leapfrog')
[149]	249	! CALL leapfrog_scheme
[129]	250	!
	251	! CASE ('adam_bashforth')
[149]	252	! CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
	253	! CALL adam_bashforth_scheme
[133]	254	CASE DEFAULT
[129]	255	STOP
	256	END SELECT
	257	END DO
[130]	258
[148]	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
[130]	263
[133]	264	IF(MOD(it+1,itau_adv)==0) THEN
[148]	265	! CALL transfert_request(f_wfluxt,req_i1) ! FIXME
[149]	266	! CALL transfert_request(f_hfluxt,req_e1) ! FIXME
[138]	267
[135]	268	CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step
[134]	269	fluxt_zero=.TRUE.
[138]	270
	271	! FIXME : check that rhodz is consistent with ps
[148]	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
[133]	285	END IF
[149]	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	!----------------------------------------------------
[124]	294	! CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
[129]	295	ENDDO
	296
[12]	297	CONTAINS
	298
[130]	299	SUBROUTINE Euler_scheme(with_dps)
[12]	300	IMPLICIT NONE
[130]	301	LOGICAL :: with_dps
	302	INTEGER :: ind
[148]	303	INTEGER :: i,j,ij,l
[145]	304	CALL trace_start("Euler_scheme")
	305
[130]	306	DO ind=1,ndomain
[138]	307	CALL swap_dimensions(ind)
	308	CALL swap_geometry(ind)
[130]	309	IF(with_dps) THEN
[148]	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))
[130]	321	END IF
[148]	322
[130]	323	u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
	324	du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
[148]	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
[130]	337	ENDDO
[133]	338
[145]	339	CALL trace_end("Euler_scheme")
	340
[12]	341	END SUBROUTINE Euler_scheme
[120]	342
[129]	343	SUBROUTINE RK_scheme(stage)
[120]	344	IMPLICIT NONE
	345	INTEGER :: ind, stage
[129]	346	REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
[120]	347	REAL(rstd) :: tau
[148]	348	INTEGER :: i,j,ij,l
[145]	349
	350	CALL trace_start("RK_scheme")
[120]	351
	352	tau = dt*coef(stage)
[129]	353
[120]	354	DO ind=1,ndomain
[138]	355	CALL swap_dimensions(ind)
	356	CALL swap_geometry(ind)
[120]	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)
[129]	360
	361	IF (stage==1) THEN ! first stage : save model state in XXm1
[148]	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
[120]	382	END IF
[129]	383	! updates are of the form x1 := x0 + tau*f(x1)
[148]	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
[133]	403	IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
	404	hflux=f_hflux(ind); hfluxt=f_hfluxt(ind)
[138]	405	wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
	406	CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
[133]	407	END IF
[120]	408	END DO
[145]	409
	410	CALL trace_end("RK_scheme")
	411
[120]	412	END SUBROUTINE RK_scheme
	413
[12]	414	SUBROUTINE leapfrog_scheme
	415	IMPLICIT NONE
	416	INTEGER :: ind
	417
[145]	418	CALL trace_start("leapfrog_scheme")
	419
[12]	420	DO ind=1,ndomain
[138]	421	CALL swap_dimensions(ind)
	422	CALL swap_geometry(ind)
[12]	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(:)+2dtdps(:)
	439	u(:,:)=um2(:,:)+2dtdu(:,:)
	440	theta_rhodz(:,:)=theta_rhodzm2(:,:)+2dtdtheta_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
[145]	447
	448	CALL trace_end("leapfrog_scheme")
	449
[12]	450	END SUBROUTINE leapfrog_scheme
	451
[129]	452	SUBROUTINE leapfrog_matsuno_scheme(stage)
[12]	453	IMPLICIT NONE
[129]	454	INTEGER :: ind, stage
	455	REAL :: tau
[145]	456
	457	CALL trace_start("leapfrog_matsuno_scheme")
	458
	459	tau = dt/nb_stage
[12]	460	DO ind=1,ndomain
[138]	461	CALL swap_dimensions(ind)
	462	CALL swap_geometry(ind)
	463
[12]	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
[129]	470	! IF (MOD(it,matsuno_period+1)==0) THEN
	471	IF (stage==1) THEN
[12]	472	psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
[129]	473	ps(:)=psm1(:)+tau*dps(:)
	474	u(:,:)=um1(:,:)+tau*du(:,:)
	475	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
[12]	476
[129]	477	! ELSE IF (MOD(it,matsuno_period+1)==1) THEN
	478	ELSE IF (stage==2) THEN
[12]	479
[129]	480	ps(:)=psm1(:)+tau*dps(:)
	481	u(:,:)=um1(:,:)+tau*du(:,:)
	482	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
[12]	483
	484	psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
	485	psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)
	486
	487	ELSE
	488
[129]	489	ps(:)=psm2(:)+2taudps(:)
	490	u(:,:)=um2(:,:)+2taudu(:,:)
	491	theta_rhodz(:,:)=theta_rhodzm2(:,:)+2taudtheta_rhodz(:,:)
[12]	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
[145]	499	CALL trace_end("leapfrog_matsuno_scheme")
[12]	500
	501	END SUBROUTINE leapfrog_matsuno_scheme
	502
	503	SUBROUTINE adam_bashforth_scheme
	504	IMPLICIT NONE
	505	INTEGER :: ind
	506
[145]	507	CALL trace_start("adam_bashforth_scheme")
	508
[12]	509	DO ind=1,ndomain
[138]	510	CALL swap_dimensions(ind)
	511	CALL swap_geometry(ind)
[12]	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(23dps(:)-16dpsm1(:)+5dpsm2(:))/12
	523	u(:,:)=u(:,:)+dt(23du(:,:)-16dum1(:,:)+5dum2(:,:))/12
	524	theta_rhodz(:,:)=theta_rhodz(:,:)+dt(23dtheta_rhodz(:,:)-16dtheta_rhodzm1(:,:)+5dtheta_rhodzm2(:,:))/12
	525
	526	dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
	527	dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
	528
[50]	529	ENDDO
[145]	530
	531	CALL trace_end("adam_bashforth_scheme")
[12]	532
	533	END SUBROUTINE adam_bashforth_scheme
	534
[50]	535	END SUBROUTINE timeloop
[133]	536
[138]	537	SUBROUTINE compute_rhodz(comp, ps, rhodz)
[133]	538	USE icosa
	539	USE disvert_mod
[138]	540	LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
[133]	541	REAL(rstd), INTENT(IN) :: ps(iim*jjm)
[146]	542	REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
[138]	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
[133]	553	DO l = 1, llm
[138]	554	DO j=jj_begin-dd,jj_end+dd
	555	DO i=ii_begin-dd,ii_end+dd
[133]	556	ij=(j-1)*iim+i
[149]	557	m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
[138]	558	IF(comp) THEN
	559	rhodz(ij,l) = m
	560	ELSE
	561	err = MAX(err,abs(m-rhodz(ij,l)))
	562	END IF
[133]	563	ENDDO
	564	ENDDO
	565	ENDDO
[138]	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
[149]	572	PRINT *, 'No discrepancy between ps and rhodz detected'
[138]	573	END IF
	574	END IF
	575
[133]	576	END SUBROUTINE compute_rhodz
	577
[138]	578	SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
[133]	579	USE icosa
	580	REAL(rstd), INTENT(IN) :: hflux(3iimjjm,llm), wflux(iim*jjm,llm+1)
	581	REAL(rstd), INTENT(INOUT) :: hfluxt(3iimjjm,llm), wfluxt(iim*jjm,llm+1)
	582	REAL(rstd), INTENT(IN) :: tau
[134]	583	LOGICAL, INTENT(INOUT) :: fluxt_zero
[148]	584	INTEGER :: l,i,j,ij
	585
[134]	586	IF(fluxt_zero) THEN
[138]	587	! PRINT *, 'Accumulating fluxes (first)'
[134]	588	fluxt_zero=.FALSE.
[148]	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
[134]	608	ELSE
[138]	609	! PRINT *, 'Accumulating fluxes (next)'
[148]	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
[134]	630	END IF
[133]	631	END SUBROUTINE accumulate_fluxes
[12]	632
	633	END MODULE timeloop_gcm_mod

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: