Context Navigation

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

Last change on this file since 159 was 159, checked in by dubos, 11 years ago
Towards Lagrangian vertical coordinate (not there yet)
File size: 18.4 KB

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

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

Download in other formats: